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SUMMARY 

Helicopters can experience high vibration levels, which reduce passenger comfort and 
cause progressive damage to the aircraft structure and on-board equipment. Because the 
primary source of excitation is typically the main rotor, special rotor control systems 
have been proposed to reduce vibrations at the source. This study addresses one such 
system, generally known as “Higher Harmonic Control” (HHC), because it consists of 
superimposing high frequency rotor inputs to the conventional low frequency ones used to 
control and maneuver the helicopter. Because both the primary flight control system and 
the HHC system act on the main rotor, the risk of adverse interactions between the two 
systems exists. This research focuses on these interactions, which have never been studied 
before due to the lack of suitable mathematical models. 

The key ingredient is an accurate linearized model of the helicopter, which includes the 
higher harmonic rotor response, and both the Automatic Flight Control System (AFCS) 
and the HHC system. Traditional linearization techniques lead to a system with periodic 
coefficients. Although Floquet theory can be used to study such periodic systems, there 
are far more control system design theories and software tools available for linear time- 
invariant systems than for periodic systems. Additionally, the theoretical evaluation of the 
helicopter handling qualities requires linear time-invariant systems. 

This research describes a new methodology for the extraction of a high-order, linear time 
invariant model, which allows the periodicity of the helicopter response to be accurately 
captured. This model provides the needed level of dynamic fidelity to permit an analysis 
and optimization of the AFCS and HHC algorithms. The key results of this study indicate 
that the closed-loop HHC system has little influence on the AFCS or on the vehicle handling 
qualities, which indicates that the AFCS does not need modification to work with the HHC 
system. On the other hand, the results show that the vibration response to maneuvers must 
be considered during the HHC design process, and this leads to much higher required HHC 
loop crossover frequencies. This research also demonstrates that the transient vibration 
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Moffett Field, California 
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responses during maneuvers can be reduced by optimizing the closed-loop higher harmonic 
control algorithm using conventional control system analyses. 
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1 Introduction 


1.1 Motivation 

Excessive vibration levels can reduce mission effectiveness on military aircraft and 
decrease passenger comfort and acceptance on commercial aircraft. Even moderate 
fuselage vibrations reduce the reliability of on board equipment, such as avionics (ref. 1). 
Maintenance costs can be significantly reduced if airframe vibrations are reduced. It 
has been estimated (ref. 2) that by reducing the fuselage vibrations in the Sikorsky 
UH-60 helicopter from 0.2g to O.lg, $80,000 per aircraft per year can be saved in 
direct maintenance costs. This is a savings of $160 million/year for a fleet of 2,000 
aircraft! These savings are achieved primarily from reduced component failures due to 
vibration. Consequently, vibration reduction is a high priority for helicopter designers and 
manufacturers. 

The major source of vibration is the unsteady aerodynamic environment experienced by 
the rotor blades including blade/vortex interaction, retreating blade stall, and blade/fuselage 
aerodynamics interaction. These blade loads are then transmitted through the hub, resulting 
in vibration of the elastic fuselage. The traditional approaches for reducing helicopter 
vibration are generally passive methods. They attack the vibration problem by increasing 
the number of blades, isolating the transmission system, applying hub absorbers, installing 
bifilars, or adding dynamic absorbers. These systems are heavy and have narrow frequency 
effectiveness ranges. Over the past two decades, the helicopter industry, government 
and academia have demonstrated that Higher Harmonic Control (HHC) is an effective 
method for vibration reduction. HHC technology may be able to achieve greater vibration 
reduction with less weight than traditional approaches by suppressing vibration at the 
source. Typically, the HHC input frequency has been nl rev, where n is the number of rotor 
blades, but other frequencies have also been utilized. A detailed survey of the extensive 
work in the area has been presented by Friedmann (ref. 3) and Teves et al. (ref. 4). 

The active rotor control system for vibration suppression is shown in figure 1.1. The 
helicopter control system generally consists of two control systems: Automatic Flight 
Control System (AFCS) and HHC system. In figure 1.1, the AFCS manages the helicopter 
stability and controls, and the HHC system suppresses the helicopter vibration. The HHC 
loop consists of three basic components. First, the data acquisition system (A/D, analog- 
to-digital converter) receives the helicopter hub loads Z (t) and converts them to the digital 
signal Z(k). Next, the harmonic analyzer extracts the nl rev harmonic components of the 
hub loads and forwards them to the HHC controller. Fast, the HHC controller computes 
the ideal HHC inputs 9{k) for vibration suppression. 

The rotor control system does not receive the new HHC input from HHC controller 
at every time step. The HHC input update rate (number of HHC input update per rotor 
revolution) depends on the time required to complete the data acquisition and post data 
processing, and has a strong influence on the HHC loop stability margin. HHC update 
rates from 0.5 to 16 times per rotor revolution have been implemented on several wind 
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Figure 1.1. Typical active rotor control system for vibration suppression. 


tunnel and flight tests (refs. 5-10); however, the once-per-revolution (or 1/rev) HHC update 
is most commonly implemented. To date, very little information is available on potential 
interactions between the HHC and AFCS. Published literature describes results from flight 
tests, wind tunnel tests, and numerical simulations with either closed-loop HHC or closed- 
loop AFCS, but not with both types of loops closed simultaneously. Because of the periodic 
nature of the helicopter, HHC is a control system application that has developed without 
the benefit of standard control system analysis techniques. Wereley and Hall (ref. 11) have 
studied the stability of the closed-loop HHC system, but the plant model was assumed 
to be quasi-static, and did include periodic behavior of the rotor system. Therefore, the 
achievable bandwidth of HHC algorithms was limited by the quasi-static assumption on the 
plant model. The HHC performance improvement could only be achieved by including the 
periodic behavior of the rotor system in the plant model and developing a control algorithm 
for the periodic time-varying plant model. Although Floquet theory can be used to study 
the periodic time-varying system, there are far more control system design methods and 
software tools available for linear time invariant systems than for periodic systems. 

Furthermore, the effects of the HHC system on vehicle handling-qualities and 
maneuverability remained unknown. There are several analyses that are important for 
evaluating the handling-qualities of the helicopter system that currently cannot be carried 
out. These include calculations of gain and phase margins with the closed-loop HHC and/or 
AFCS, crossover frequency of the HHC loops, and closed-loop stability of helicopter 
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dynamics with closed-loop HHC. These quantities can be easily obtained from a linear 
time-invariant system. Therefore, there is a need for linear time-invariant approximations 
that can accurately model the coupled rotor-fuselage dynamics, including the higher 
harmonic response of the rotor. Such time-invariant linearized approximation methods 
are not currently available. 

1.2 Literature review 

1.2.1 Higher harmonic control technology 

In 1952, Stewart (ref. 12) showed the potential effectiveness of HHC in alleviating 
retreating blade stall. The use of a second harmonic control (2/rev) was shown to 
redistribute rotor disk loads and suppress retreating blade stall. By delaying the retreating 
blade stall to a higher forward speed, the speed limitation of a helicopter could be 
raised. Based on his analysis, the advance ratio could be increased by 0.1. However, 
the analysis was based on a rigid flapping blade and airloads were calculated with quasi- 
static aerodynamics and uniform inflow distribution. The transonic effects, separated flow 
conditions, unsteady aerodynamics, blade flexibility, and non-uniform inflow distribution 
were all neglected. 

In 1961, Arcidiacono (ref. 13) extended Stewart’s research by including both 2/rev and 
higher harmonic blade pitch control. The analyses showed that a combination of 2 and 
3/rev HHC inputs could be used to delay retreating blade stall to an even higher advance 
ratio than that reported by Stewart. Neither Stewart nor Arcidiacono considered the effects 
of HHC input on vibratory hub loads. 

In 1961, the first HHC flight test was carried out to investigate the feasibility of using 
HHC for vibration suppression on UH-1A 2-bladed helicopter (ref. 14). A series of flight 
tests was conducted by Bell Helicopter Company to determine the effects of HHC on rotor 
performance, blade airloads, blade loads, control loads, hub loads, and airframe vibration. 
The investigators noted that no reduction in shaft torque was observed. The investigation 
also showed that drag reduction in the retreating side of the rotor was accompanied by an 
increase in profile drag in the fore and aft portion of the rotor disk. These results confirmed 
Stewart’s finding, and indicated that 2/rev HHC input could be used to change the rotor 
disk loading. 

In 1972, McCloud (refs. 15, 16) reported the first full-scale wind tunnel investigation on 
HHC. The rotor model was a two-bladed teetering rotor with propulsive jet flaps. A large jet 
flow was expelled from the blade trailing edge to propel the rotor and the HHC was applied 
through the angular deflection of the jet flow. The experiments showed that the vibratory 
hub load reduction was accompanied by an increase in the blade bending moments. The 
HHC inputs required for the vibration suppression were found to vary with rotor forward 
speed. 

In 1975, McHugh and Shaw (refs. 17, 18) conducted a series of wind tunnel experiments 
on a four-bladed soft-inplane hingeless rotor model. The HHC was implemented in the 
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non-rotating frame, and HHC inputs were applied by oscillating the swashplate with servo- 
actuators. Results from the experiments indicated that the vibratory hub moments could 
be suppressed effectively without a significant increase in blade stresses. The experiments 
also indicated that all five components of the 4/rev hub loads (lateral, longitudinal, and 
vertical forces; pitching and rolling moments) could be reduced simultaneously with three 
HHC inputs. 

In 1979, a wind tunnel investigation of HHC on a four-bladed hingeless rotor model was 
conducted by Shaw and Albion (refs. 19,20) in the Boeing V/STOL Wind Tunnel. The 
rotor model was Mach scaled and operated at full-scale tip speed. The HHC inputs were 
applied through swashplate excitation. The closed-loop HHC controller simultaneously 
reduced the 4/rev vertical hub shear, hub pitching and rolling moments by up to 90%. 
The closed-loop transient behaviors were studied by introducing a step disturbance in the 
swashplate command. The results showed that the disturbance was suppressed within two 
rotor revolutions, which confirmed the quasi-static assumption made on the HHC model. 

In 1980, Shaw (ref. 21) presented the results of a comprehensive analytical investigation 
of HHC. He compared the potential benefits of servo-flap versus conventional blade root 
feathering and studied the automatic in-flight adaptive algorithm. The investigation was 
based on a coupled modal analysis and included a vortex wake induced flow calculation. 
An approximation to the Theodorsen lift deficiency function was used to include the effect 
of the shed wake. A transfer matrix (T-matrix) approach was implemented to relate 
the higher harmonic hub loads to the HHC inputs. The analytical results showed that 
nonlinearities in the HHC input-output model were small. The results also indicated that the 
vibration suppression was caused by mutual cancellation between aerodynamic and inertial 
components of the transmitted vibratory loads at the blade roots. With the HHC inputs the 
control loads were increased by roughly 30%, and the change in rotor performance was 
negligible. For changing flight conditions, the closed-loop HHC controller with fixed gain 
performed satisfactorily over an advance ratio range of 0.2. An adaptive gain controller 
was used in cases where the fixed gain controller performed poorly. For the adaptive gain 
controller, the model parameters were estimated using a Kalman filter. Simulation results 
showed that the adaptive controller performed well for varying flight conditions. 

In 1981, Molusis, Hammond and Cline (ref. 22) studied several HHC algorithms 
for vibration suppression, and the algorithm performance was evaluated in wind tunnel 
testing. The rotor model was a Mach-scaled four-bladed articulated rotor. The HHC 
controllers were configured to suppress the 4/rev vertical, longitudinal, and lateral signals 
from a triaxial accelerometer mounted beneath the rotor in the non-rotating frame. The 
advance ratio was varied between 0.2 to 0.4. The HHC system was modeled using a 
T-matrix approach. The HHC algorithms studied were separated into two groups: the 
adaptive controllers and the gain-scheduled controllers. Each type of controller was 
further classified. The adaptive controllers were classified into deterministic controllers 
and cautious controllers. The gain-scheduled controllers were classified into perturbation 
controllers and proportional controllers. The wind tunnel results showed that the gain- 
scheduled controllers performed poorly, possibly due to the nonlinear behavior of the 
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HHC model. The deterministic (adaptive) controller was shown to significantly reduce 
the steady-state vibration levels, but there were large transient responses that occurred 
before the vibration converged to the steady-state value. The authors noted that the cautious 
controller offered the best performance among the four controllers. 

In 1981, the performance of four different feedback controllers or regulators were 
investigated by Chopra and McCloud (ref. 8) for the multicyclic control of helicopter 
vibration. These controllers were open-loop and closed-loop with off-line and on-line 
identification. The off-line identification of model characteristics was made using the least- 
squared-error method and used a succession of input and output measurements. The on-line 
identification of model characteristics was computed using a Kalman-filter solution. The 
optimal controls were calculated by minimizing the quadratic performance function based 
on response and control inputs. Both global (linear) and local (piecewise linear) models 
were simulated. The results showed that the closed-loop controller with a local model 
using on-line identification techniques performed the best. For the cases with large initial 
errors in the transfer matrix, large overshoots were found in the transient response using 
this controller. 

In 1982, Wood et al. (refs. 23, 24) conducted a HHC flight test on a modified Hughes 
OH-6 A helicopter with a gross weight of about 3,000 lb. The HHC input was achieved by 
blade root feathering through the 4/rev swashplate oscillations. A triaxial accelerometer 
was mounted beneath the pilot seat to sense and feed back the 4/rev vibrations to the 
HHC controller. The aircraft was flown from hover to 100 knots with the HHC system 
operated in open-loop (manually) and closed-loop (computer controlled). For the closed- 
loop controller, the cautious controller presented in reference 22 was used. The test results 
indicated that up to 90% reduction in vibration could be obtained with HHC amplitude less 
than 1°. 

During the 1980s, extensive research on the use of HHC implemented in the form of 
Individual Blade Control (IBC) was carried out by Ham and his coworkers (refs. 25,26). 
The potential applications of IBC that were proposed included reducing the severe effects 
of atmospheric turbulence, retreating blade stall, blade-vortex interaction, blade-fuselage 
interference, and blade instabilities, while providing improved flighting qualities and 
automatic blade tracking. The theoretical analysis showed that the rotor blade flapping, 
inplane, and torsional motion could be reduced by feedback control of the effective inertia, 
damping, and stiffness of the appropriate modes. 

In 1985, Shaw et al. (ref. 6) described wind tunnel tests on a dynamically scaled three- 
bladed CH-47D Chinook rotor in the Boeing V/STOL Wind Tunnel. The 2, 3, and 4/rev 
HHC inputs were applied to suppress the 3/rev vertical hub force and the 2 and 4/rev 
rotating inplane hub shears throughout a wide test envelope which included trimmed flight 
up to 188 knots. The open-loop tests were conducted to obtain transfer matrices under 
several flight conditions. These transfer matrices were used with fixed- or gain-scheduled 
controllers. The wind tunnel results showed that a fixed-gain controller with a local model 
can suppress more than 90% in all three vibratory hub shear components. The wind tunnel 
results indicated that the gain-scheduled controller performed as well as the fixed-gain 
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controller. The adaptive controllers, similar to those in reference 22, were either unstable 
or ineffective in suppressing the vibratory loads. 

In 1986, Polychroniadis and Achache (refs. 27,28) discussed the application of HHC on 
an Aerospatiale SA-349 Gazelle helicopter (4,500 lb) for vibration reduction and noise 
reduction, and included a performance analysis based on both theoretical studies and 
wind tunnel testing. The HHC input was achieved by blade root feathering through the 
4/rev swashplate oscillations. The HHC controller was a self-adaptive controller that used 
vibration sensors placed at pre-selected locations in the aircraft. The test results showed a 
70 to 90% reduction in vibration was achieved at forward speeds up to 135 knots. 

In 1994 and 1995, Jacklin et al. (refs. 29, 30) described two wind tunnel tests that 
evaluated the effects of IBC at various frequencies on rotor performance, vibrations, and 
acoustics using a full-scaled BO-105 helicopter rotor. The IBC system, developed by ZF 
Luftfahrttechnik, was tested on the NASA/ Army Rotor Test Apparatus in the NASA Ames 
40- by 80-Foot Wind Tunnel. Test results indicated that a single-frequency IBC input of 
2-4/rev could simultaneously reduce all 4/rev rotor balance forces and moments by up to 
70% at 43 knots. 

In 2002, U. T. P. Arnold (ref. 10) described the certification, ground and flight testing of 
an experimental IBC system for a Sikorsky CH-53G helicopter with a gross weight of about 
68,000 lbs. The primary goal of the IBC system was to extend the service life of the CH- 
53 by reducing the component fatigue and failure induced by high vibratory stresses. The 
IBC system was designed, built, installed, and certified by ZF Luftfahrttechnik, GmbH. The 
IBC system, weighing less than 1% of the helicopter maximum take-off weight, completely 
integrated all mechanical and hydraulic components into the rotating frame. The IBC 
controller was based on a second order T-matrix model. Initial test results showed a high 
effectiveness of IBC in reducing vibration with a relatively small single harmonic input of 
±0.15°. 

Most of the active vibration control algorithms discussed above were implemented in 
frequency domain. In 1980, Du Val and Gupta (ref. 31) proposed a time domain approach 
for the active control of helicopter vibration. The controller was designed by optimizing a 
cost function, which placed a large penalty on fuselage vibration at n/rev frequency. The 
fuselage accelerations were passed through an undamped second-order system tuned to the 
n/rev frequency. At the resonant frequency, the regulator locked onto the magnitude and 
phase of the fuselage accelerations without the need for harmonic analysis. By placing 
an infinite weighting on the n/rev response, a controller is guaranteed to drive the n/rev 
response to zero. Because the dynamics of the rotor and fuselage were included in the 
plant model, the quasi-static assumption was no longer necessary. Also, because the state 
feedback was used, on-line identification of the model parameters was not necessary. This 
method assumed the system was linear time-invariant (LTI), not linear time-periodic (LTP). 
The standard linear analysis techniques and software tools could therefore be applied. 

In 1989, Wereley and Hall (refs. 1 1, 32) presented a framework to provide the evaluation 
of active vibration control algorithms performance in terms of classical control theory. 
They showed that HHC was fundamentally similar to the sinusoidal disturbance rejection 
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techniques of classical control. By treating the periodic disturbance as a stochastic rather 
than a deterministic phenomenon, the methods of Shaw et al. (ref. 6) and Du Val and 
Gupta (ref. 31) could be compared quantitatively. The authors indicated that the achievable 
bandwidth of HHC algorithms was limited by the quasi-static and linear time-invariant 
assumptions on the plant model. The HHC performance improvement could only be 
achieved by including the periodic behavior of the rotor system in the plant model and 
developing a control algorithm for the periodic time- varying plant model. 

In 2000, Spencer (ref. 33) presented the open and closed-loop wind tunnel testing of a 
Mach-scaled active rotor system with piezoelectric bender actuated trailing-edge flaps for 
active vibration control. The closed-loop vibration suppression tests were conducted at 
several advance ratios and collective settings. The controller design is based on a radial 
basis neural network which is used to approximate the command input to the active rotor. 
The controller is implemented in discrete time by sampling the hub loads and control input 
at 1/rev frequency. The optimum set of network weights is determined by minimizing 
the cost function which is based on the vibration response and command input. One of 
the advantages of the neural network controller is that it simultaneously learns while it 
commands the on-blade actuator, thus adaptively suppressing the blade or hub vibrations. 
No off-line training of the network is required. These tests successfully reduced the 
4/rev oscillatory fixed frame thrust, pitching moment, and rolling moment levels up to 
90%. A transient vibration control test was also conducted by varying the rotor speed, 
wind speed, and the collective pitch angle to simulate maneuvering flight. For all three 
individual perturbations, the neural-controller was unable to compensate vibration response 
fluctuation. The authors indicated the controller was not able to react fast enough to the 
perturbations because of hardware limitations. 

1.2.2 Linear models 

In 1981, Howlett (ref. 34) presented a nonlinear mathematical model known as GenHel, 
based on the Sikorsky UH-60A Black Hawk helicopter, for performance and handling- 
qualities evaluations. The rotor was modeled with a rigid blade flap and rigid blade lag 
degree of freedom. The torsional dynamics were modeled using a simple dynamic twist 
model. The aerodynamic forces on the rotor were computed using blade element theory 
and quasi-static aerodynamics. Aerodynamics coefficients of the blade were provided by 
the look-up tables as a function of the angle of attack and Mach number. The fuselage 
was modeled as a rigid body with aerodynamic coefficients of the fuselage and empennage 
provided by the look-up tables as a function of angle of attack. 

The GenHel simulation model could provide a linear model, but it was limited to six 
fuselage degrees of freedom. The linearization was performed numerically by perturbing 
each of the states and controls, and using finite difference approximations. Because of 
the unusual flight dynamic model implementation in GenHel, the perturbation scheme was 
not straightforward. For instance, the fuselage states and control inputs were perturbed 
one at a time about the trim condition to produce a 9-state rigid body linear model from 
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the nonlinear GenHel mathematical model. The rotor equations of motion continued to 
be integrated while all rigid body acceleration integrations were suppressed. The change 
in the state derivatives was calculated when the rotor response had reached a steady-state 
condition. This method produced a six-fuselage-degree-of-freedom linear model with a 
quasi-static rotor; i.e., the dynamics of the rotor system were not modeled. 

In 1982, Zhao and Curtiss (ref. 35) developed the first linear model of a helicopter that 
included blade dynamics for forward flight. This linear model had 24 or 27 states depending 
on whether dynamic inflow was included. Flap and lag degrees of freedom were modeled 
by transforming the rotor equations of motion into the non-rotating frame using multi-blade 
coordinate transformation. Only the collective and first two cyclic modes for each rotor 
degree of freedom were retained. Unsteady aerodynamic effects were introduced through 
the dynamic inflow model. A flat vortex wake model was used to approximate the effects 
of the main rotor wake interference on the tail surfaces and tail rotor. The linear model 
was derived using a symbolic mathematic manipulation program to obtain an analytical 
solution. The linearization process consisted of expressing the time dependent variables in 
the equations of motion as the sum of the trim value and time dependent perturbation about 
the trim value. A linear model could be obtained by applying order reduction and setting 
the remaining perturbation quantities to zero. 

In 1986, Chen and Tischler (ref. 36) discussed the method of developing the simplified 
analytical linearized model from the flight test data by using modem system or parameter 
identification techniques. The simplified analytical model could be used for handling- 
qualities evaluation, design of stability and control augmentation systems, and ground 
simulator validation. Authors stated that the importance of recognizing that each lower- 
order model used for rotorcraft parameter identification had a limited range of applicability. 
They also discussed the benefits and limitations of using frequency sweeps as flight test 
input signals for identification of frequency response for rotorcraft and for the subsequent 
fitting of parametric transfer function models. The authors concluded that analytical 
modeling and understanding the limitation of lower-order models could be more important 
than merely relying on the identification algorithms. 

In 1987, Miller and White (ref. 37) presented an algorithm called Exponential Basis 
Function (EBF) which allowed computer generation of a comprehensive coupled rotor- 
fuselage nonlinear model. EBF represented the position vector of a generic mass element 
of the helicopter exponentially, and was used to simplify the differentiation of the position 
vector. EBF was used to write the time dependent coordinate transformation as the product 
of constant matrices and matrix exponentials. Since the multiplication of exponentials 
is equivalent to addition of exponential arguments, multiplication of the transformation 
matrices could be accomplished by adding matrix exponentials. The transformation 
matrices written in EBF could also be differentiated easily. The equations of motion were 
formulated using Lagrange’s equation. The rotor degrees of freedom were transformed to 
a non-rotating frame using the multi-blade coordinate transformation. The rotor dynamics 
included rigid body flap and lag degrees of freedom. Engine rotor speed, fuselage rigid 
body degrees of freedom, and the inflow dynamics were also modeled. Authors stated 
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that the linear model obtained through EBF could be used to analyze handling-qualities 
phenomena for highly augmented helicopters when realistic trim conditions and high-order 
dynamics were considered. 

In 1990, Kim, Celi, and Tischler (ref. 38) developed a high-order linearized model 
of helicopter flight dynamics extracted from a nonlinear time domain simulation. The 
model had 29 states that described the fuselage rigid body degrees of freedom, the flap 
and lag dynamics in a non-rotating frame, the inflow dynamics, and the delayed entry of 
the horizontal tail into the main rotor wake. The blade torsional degree of freedom was 
approximated using a pseudo-modal approach. In GenHel, the calculation of forces and 
moments acting on the helicopter at a given instant in time was solved sequentially; the 
rotor equations of motion were solved first, and the fuselage equations of motion were 
solved next. Because of this separation, the equations of motion were not perturbed 
simultaneously, which could cause inaccuracies in the solution at higher frequencies. The 
perturbation process was also complicated by this splitting solution process. Therefore, 
GenHel could only produce a six-fuselage-degree-of-freedom linear model with a quasi- 
static rotor. 

To carry out a theoretically rigorous linearization and retain the rotor dynamics within 
the linear model, the mathematical model of the helicopter as implemented in GenHel 
was extensively modified to a first-order, state variable form. This required several 
modifications including solving both the rotor and the fuselage equations of motion 
simultaneously. The linear model was validated against the nonlinear model, and the results 
showed a good agreement between these two models for small amplitude control inputs. In 
case of large amplitude inputs, which violated the small perturbation assumption inherently 
contained in the linear model, the agreement deteriorated greatly. 

1.3 Objectives of study 

The objectives of this study are as follows: 

• Develop a methodology for the derivation of linearized, time-invariant, state- 
space models of coupled rotor-fuselage dynamics that include the effects of higher 
harmonic response of rotor and fuselage to both higher harmonic pitch control and 
pilot inputs. 

• Apply the new linear state-space models for a study of the potential interactions 
between a higher harmonic control system and an automatic flight control system, 
including any impact on handling-qualities. 

It should be pointed out that this research does not focus on the method to improve 
the helicopter vibratory hub load predictions. A comprehensive analysis on this topic 
is beyond the scope of the present study. The helicopter simulation model used in this 
study is adequate to capture the first-order effects, but it may not be sufficient for accurate 
quantitative predictions of vibratory hub loads. 
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1.4 Principal contributions 

• Implemented a HHC in a flight dynamics model for a free flight condition to 
investigate the interaction between HHC and AFCS. 

• Developed a linear time-invariant state-space approximation that accurately models 
the coupled rotor-fuselage dynamics including the higher harmonic response of the 
rotor. This coupled high-order linear model provides the needed level of dynamic 
fidelity to permit study of AFCS and HHC interaction. 

• Provided detailed analyses on the HHC/AFCS interaction, and developed an 
improved HHC controller to reduce the vibration transients during the maneuvering 
flight. 

1.5 Organization of the document 

Chapter 2 describes the mathematical model of the helicopter and provides the solution 
method for the trim calculation, linearization, time integration, and vibratory hub 
load calculation. 

Chapter 3 is devoted to the HHC system for the vibration suppression. The inner working 
of the harmonic analyzer, HHC controller, and the HHC update scheme are discussed 
in detail. The methods of obtaining the continuous-time domain equivalent for each 
component are also presented. 

Chapter 4 presents a new linearization method that converts a nonlinear system to a linear 
time-invariant system while capturing the n/rev characteristic of the helicopter. The 
new linear model was validated by comparing vibratory hub loads and rotor states 
for both higher harmonic inputs and piloted input at several forward speeds. 

Chapter 5 presents the results of the HHC/AFCS interaction study. The effect of HHC 
input on handling-qualities was tested for both open-loop and closed-loop HHC 
systems. This chapter also discusses the effect of the HHC on the vibration transients 
during maneuvers, and develops a new HHC algorithm to overcome the problem. 

Chapter 6 presents conclusions of the study and recommendations for future work. 
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2 Mathematical Model 


This chapter contains a brief history of the helicopter mathematical model used in this 
study, followed by the definition and implementation method of the HHC system, the 
methods used to calculate the helicopter trim states, a linearized model, and the time 
history. Next, the definition and the computation method of the vibratory hub load are 
discussed. Finally, the last section presents the method of extracting n/rev vibration using 
Fourier series approximation. 

2.1 History of helicopter simulation model 

The flight dynamic simulation model used in this study is originally from the helicopter 
simulation model GenHel (ref. 34) specialized for the Sikorsky UH-60 Black Hawk. The 
rotor was modeled with a rigid blade flap and rigid blade lag degrees of freedom. The 
torsional dynamics were modeled using a pseudo-modal approach. The fuselage was 
modeled as a rigid body with aerodynamic coefficients of the fuselage and empennage 
provided by the look-up table. The fidelity of GenHel model was improved by 
B allin (ref. 39) who also implemented the engine model. Kim (refs. 40,41) included 
the main rotor inflow model using the Pitt-Peters dynamic inflow model (ref. 42). A 
new trim procedure was also developed with the equations of motion presented in first- 
order state-space form. This allows a linear time-invariant model to be extracted using a 
perturbation-averaging method. The model developed by Kim was named UM-GenHel. 
UM-GenHel was continued in a series of calibrations based on the flight test data at NASA 
Ames Research Center. This version of UM-GenHel was renamed FORECAST, and is 
widely used in flight dynamics analysis at NASA Ames Research Center. 

At the same time, the UM-GenHel remained at the University of Maryland as a research 
helicopter model. Tumour (ref. 43) extended the rotor blade modeling in UH-GenHel 
by including the aeroelastic rotor, which was originally developed by Celi (ref. 44) and 
extended by Spence (ref. 45) to include the coupled rotor/fuselage formulation. Tumour 
also added the finite state wake (ref. 46) and the Leishman-Nguyen (ref. 47) state-space 
unsteady aerodynamics model. This version of the research model was renamed by Tumour 
as FlexUM. Theodore (ref. 48) extended the inflow flow model to include the maneuvering 
Free Wake model (ref. 49), which improves the off-axis response predictions. A full BO- 
105 helicopter configuration is also added to the FlexUM. The research model was renamed 
to HeliUM by Theodore. 

2.2 Helicopter model 

The basic formulation and solution of the equations are unchanged with respect to the 
previous works. The helicopter model used in this study is similar to the Sikorsky UH- 
60 Black Hawk with the following simplifying assumption. The helicopter equations of 
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motion are based on a set of coupled nonlinear rotor-fuselage equations in first-order, state- 
space form. The rigid body dynamics of the helicopter is modeled using non-linear Euler 
equations. The aerodynamic coefficients of fuselage and tail surfaces are provided in the 
form of look-up tables. The blade is assumed to be straight, i.e., with zero tip sweep. The 
blade dynamics consists of two rigid blade degrees of freedom plus first torsional degree 
of freedom. The aerodynamic coefficients of the blade are also provided in the form of 
look-up tables as a function of angle of attack and Mach number. Unless stated otherwise, 
the main rotor inflow is calculated using a three-state dynamic inflow model, which yields 
linear inflow distributions over the rotor disk. Tip losses are taken into account by assuming 
that the outboard 3% of the blade does not generate lift. A one-state dynamic inflow model 
is used for the tail rotor. Stall and compressibility effects are incorporated in a quasi-static 
form, and unsteady aerodynamic effects have been neglected. Two additional assumptions 
are that the rotor speed is constant and that there is no limitation on the power supplied 
by the engine. All the results presented in this study are obtained from a coupled rotor- 
fuselage trim procedure simulating free flight conditions. All trim calculations include the 
HHC input, if one is present. In all the parametric studies, the helicopter is retrimmed every 
time the magnitude or phase of the nl rev input changes. 

2.3 HHC implementation 

The higher harmonic control inputs are implemented by varying the blade pitch at blade 
root. Unlike the real active pitch links system, stiffness of the pitch link is assumed to be 
infinitely stiff and dynamics of the active pitch links is ignored. The geometric pitch angle 
6q of the blade is given by: 

OgW = #0 + 9 lc cos O + Asp) + 6 ls sin(^ + A S p) + 6> n (^) (2.1) 

where 9 0 , 9 lc , and 9 ]s are respectively the collective, lateral cyclic, and longitudinal cyclic 
pitch, A S p is the swashplate phasing angle, A SP = —9.7°, and 9 n {^) is the 77 ,/rev input, 
defined as: 

9 n (ip) = A n cos (nip + 4> n ) (2.2) 

where A n and o n are the magnitude and phase of the n/rev input. 

2.4 Solution methods: trim 

This section presents methods to calculate the helicopter trim states. The flight condition is 
assumed to be a steady coordinated helical turn. Straight level flight then becomes a special 
case where both flight path angle and turn rate equal zero. The helicopter trim equations 
were originally developed by Chen (ref. 50), and later extended by Celi (ref. 51) to include 
the steady state response of the rotor. They are modified further by Kim (ref. 40) to consider 
the periodicity of both rotor and fuselage motion. The trim states are generally obtained 
from an algebraic trim procedure. 
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2.4.1 Algebraic trim 

The typical trim solution is based on algebraic trim. The solution of the steady state 
condition is determined by converting a set of coupled ordinary differential equations to 
a set of coupled nonlinear algebraic equations. The periodicity of the helicopter response 
must be satisfied in a steady state condition. This set of algebraic equations is then solved 
using the Powell Hybrid algorithm. The trim solution is reached when the sum of the forces 
and moments at the vehicle center of gravity are zero in one rotor revolution. 

Although this method can obtain a trim solution quickly, it does not guarantee that the 
rotor blades return to the same position after one revolution of time integration. In other 
words, time integration starting from an algebraic trim solution without control perturbation 
may not respond precisely to n-multiple/rev. This does not appear to be crucial for flight 
dynamics analyses, but has a large effect on vibration related computations. The periodic 
trim procedures can fulfill this task. 

2.4.2 Periodic trim 

There are two methods to achieve a periodic trim solution: the shooting method, and the 
time marching method. 

2.4.2.1 The shooting method 

After algebraic trim is achieved, the state vector and control vector are adjusted such that 
the state vector remains the same after integration of one rotor revolution. This is a two- 
point nonlinear boundary value problem, and is based on a shooting method (ref. 52). The 
basic idea behind the shooting method is to convert a boundary value problem (BVP) into 
an initial value problem (IVP). Given an initial guess for the parameters, an iterative solver 
is used to find values of the parameters that produce solutions that satisfy the boundary 
conditions. The method will guarantee a n/rev periodic trim, but its convergence proved 
erratic, and at least one order of magnitude more expensive computationally, compared 
with the algebraic trim procedure. 

2.4.2.2 The time marching method 

The second method is the time marching solution which is also the one used in this research. 
As stated earlier, the free flight response from the time integration starting from the 
algebraic trim solution may not have precise n-multiple/rev response. Because of unstable 
Phugoid mode, the helicopter will slowly drift away from trim. The low gain stabilization 
loop was added to ensure the helicopter does not become unstable as integration time 
increases. As the time integration continues, the n-multiple/rev response will emerge. 
Generally, the periodic trim solution can be reached within four rotor revolutions starting 
from an algebraic trim solution. At the end of time integration, the trimmed state vectors 
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around the rotor azimuth are available. Additional information about the time integration 
is presented in section 2.6. 


2.5 Solution methods: linearization of the equation of motion 

This technique consists of perturbing each state and control about an equilibrium position. 
Using this approach, the individual blade pitch is introduced in terms of the harmonics 
in the rotating frame. This method leads to systems of rotor equations containing 
periodic coefficients, which are represented in the rotating frame. The transformation 
from the rotating frame to the fixed frame is accomplished using a Multi-blade Coordinate 
Transformation (MCT, Appendix C). To remove the time dependency, the linearized 
models are computed over one rotor revolution and then averaged to obtain a LTI system in 
the fixed frame. As a consequence, this averaging eliminates the periodicity of the system 
and all the higher harmonics in both the controls and rotor response. Additional information 
about this technique is discussed in chapter 4. 


2.6 Solution methods: time integration 

The free flight response of the helicopter is computed by integrating the equations of motion 
based on a given set of initial condition and control inputs. The equations of motion are 
represented by a system of coupled nonlinear ordinary differential equations expressed 
symbolically in the first-order ODE form 

y = f(y,u;t) (2.3) 

where y is the state vector and u is the control vector. Equation 2.3 can be solved 
numerically using Adams-Bashforth method, which is a variable-step, variable-order, 
predictor-corrector, numerical method for solving linear first-order ordinary differential 
equations. It estimates the behavior of the solution curve by evaluating the derivative 
function at the old solution values along with the current solution and derivative function 
and uses the interpolation method to estimate the new solution. In other words, Adams- 
Bashforth methods try to squeeze information out of old solution points. For problems 
where the solution is smooth, these methods can be highly accurate and efficient. 

In this study, the simulation is started from the trim condition, and the equations of 
motion are integrated with respect to time. This produces time histories of all state variables 
for prescribed control inputs. Generally, control inputs include the time history of pilot 
inputs or the swashplate controls. For the HHC system, the control inputs are extended to 
the blade root pitch angle which can be prescribed as single or multiple harmonics in terms 
of the n ! rev amplitude A n and phase angle o n as stated in equation 2.2. 
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2.7 Vibration calculation 


2.7.1 Hub loads calculation 


The helicopter equations of motion can be expressed as equation 2.3. Because the coupled 
rotor-fuselage-inflow equations of motion have the state derivatives appearing on the right 
hand side of the equations, these equations are expressed as: 


y c = f c (y,y,u;t) 


(2.4) 


where y c is a vector which contains all the state derivatives appearing on the right hand 
side of equation 2.4. For instance, the flap equation for \ th blade of a simple rotor model 
(rigid flap and lag modes only) is: 

S h 


/3i = 


cos /3i{w + e [20 (p cos ^ — ? sin ^j) + p sin i/y + qcosipi] j 


+ sin Pi cos sin ^ — -ucos-^i— e (r — fl) 2 ^ 

+ cos 2 /3j| cos Q p sin ?/y + geos?/)* — 2{C t + Ft) (q sini/g 
— 2Ft sin (p sin i/ji + q cos tpi) | 

+ cos Pi sin Pi 1 2Q (r — O) - (r - Q) 2 -— Cf } 
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p COS 'Ipi 


(2.5) 


where S b and I b are the first and second blade moments of inertia about its hinge, M Aern a is 
the flap aerodynamic moment, and M LD p is the flap moments due to the lag damper. Since 
the state derivatives it, v, w, p, and q on the right hand side of equation 2.5 do not couple 
with other state derivatives, it can be rewritten as: 


where 


e = 


ih = [e] y c + (y, u; t) 


— sin / 3 cos ( cos ifj 


(2.6) 
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(2.7) 


y c = [u, v, w,p, q,r,Fl, /?i, /? 2 , /? 3 , /?4, Ci : C 2 , 03, 040 


(2.8) 


Similar expression can also be rewritten for the remaining equations (u, v, w, p, q, r, . . .). 
The resultant row vectors e are assembled into a coupling matrix E, and equation 2.4 can 
be rewritten as follow: 

y c = [E] y c + F k (y, u; t) (2.9) 
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(2.10) 


By re-arranging equation 2.9 into first-order form, y c can be solved as 

y c = [J - E] _i F k (y, u; t) 

Re-write equation 2.9 again and parse it as follow: 


y fus 


E n 

E\2 


y fus 

y mr 

c 

E 21 

E 22 


y mr 


Yc [E] y c 


+ Ftr + Ffus 

-/ 

v- 

A(y,u;t) 


y fus 

y mr 


[u : v,w : p,q,r] T 

- - 5 Pi ' P2 1 Pz i Pi : Cl t C2 * Cs • C4 


(2.11) 


(2.12) 

(2.13) 


where [Euyf us ] is the inertial acceleration due to the fuselage acceleration, [E\ 2 y mr \ is the 
inertial acceleration due to the main rotor acceleration, F mr is acceleration contributed by 
the main rotor excluding the inertial coupling term, F tr is the acceleration contributed by 
the tail rotor, and Fj vs is the acceleration contributed by the fuselage, the horizontal, and 
the vertical surfaces. Therefore, the vibratory hub loads are the sum of all the loads that are 
transmitted from the main rotor to the hub in the fixed system; i.e., F mr + |E l2 y mr ]. 


2.7.2 Cockpit vibration calculation with the rigid fuselage 

In flight test, the helicopter vibration level is measured by mounting accelerometers at 
several key areas inside the helicopter. One of the key vibration areas is the pilot station. 
The flight dynamics model (HeliUM) used in this study needed to produce the same pilot 
station acceleration in order to compare the results with the flight test data. However, 
this information is not directly available. Although HeliUM is based on a coupled 
rotor/fuselage formulation, the fuselage is actually modeled as a rigid body and does not 
contain any dynamics. All results from free flight trim procedure are only available at 
the center of gravity (CG) of helicopter. Nevertheless, the pilot station acceleration can 
be obtained by a simple transformation. Velocity at the pilot station can be expressed as 
follows: 


VpUot 

W cg + LO X R 

(2.14) 

v cg = 

[u,v,w] T 

(2.15) 

cu = 

[p, q, r} T 

(2.16) 

R = 

[x,y,z) T 

(2.17) 


where v cg is the velocity vector at CG, cu is the rotational vector of the helicopter at the 
CG, and R is the position vector from CG to the pilot station. The pilot station acceleration 
can then be calculated by differentiating equation 2.14 with respect to time, 

VpUot = v cg IwxR + wxR (2.18) 
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As the fuselage is a rigid body, the distance between the CG and the pilot station is constant; 
i.e., R = 0. Expanding the cross product term, equation 2.18 becomes 


V pilot = V C g + V X R = 


u — ry + qz 
v + r y — pz 
w — qx + py 


(2.19) 


The 4/rev vibration can then be determined by collecting v p u ot over one rotor revolution, 
and extracting its 4/rev components using the Fourier approximation. 

Figures 2.1a and 2.1b compare pilot and copilot station vibrations with flight test from 
hover to 140 kts. The flight test data represents several sets of baseline data collected over 
the span of the flight test program. The scatter in the data could be caused by changes in 
the aircraft configuration and non-ideal flight conditions during the test. 

The main rotor inflow used in HeliUM is based on a linear inflow model. The blade 
dynamics consist of a rigid blade flap, a rigid blade lag, and first blade torsion mode. The 
figures also show Yang’s (ref. 53) results from UMARC 1 , which uses 8 blade modes, a 
free wake model, and a flexible fuselage model. These two figures indicate that the cockpit 
acceleration computed from HeliUM is underestimated throughout the entire speed range. 
This study has also included additional blade flexibility (result not shown here), but the 
vibration level is very similar with the rigid blade model. 

It was believed that this under prediction could be caused by a lack of aerodynamic 
interaction. Because the linear inflow model only contains 1/rev harmonics, the higher 
harmonic airload was not excited. The linear inflow model was replaced with a free wake 
model, and the results are shown in figures 2.2a and 2.2b. As expected, the cockpit vibration 
level was greatly improved. However, the vibration level in the higher speed range is on 
the low side; especially at 120 kts, which is the baseline configuration of this research. In 
addition, the helicopter simulation with the free wake model is computationally expensive. 
The computation time required is generally over one order of magnitude higher than the 
one with a linear inflow model. 


2.7.3 Cockpit vibration calculation with the flexible fuselage 

To determine the importance of the fuselage flexibility on cockpit vibration calculations, 
the effect of the flexible fuselage is added to HeliUM. This is achieved by feeding the 
hub loads from a trim condition into a separate fully elastic 3-D fuselage model. This 
fuselage model is built using NASTRAN (ref. 54) based on a Sikorsky SH-60B (a variant 
of UH-60) helicopter fuselage. It consists of structural elements such as scalar springs, 
rods, bars, shear panels, and triangular and quadrilateral membranes for more than 8,400 
elements. NASTRAN is used to calculate fuselage mode shapes, modal mass, and stiffness. 
The resulting data are used to build a transformation matrix N which maps the 4/rev hub 
shears and moments to the cockpit accelerations at the 4/rev frequency. For example, 


'Both the flight test data and the UMARC results are the courtesy of M. Yang and I. Chopra. 


19 


the first column vector N x is obtained by applying a unit 4/rev longitudinal force F XiP 
at the NASTRAN model’s hub node and measuring all six 4/rev accelerations at the 
cockpit station. Note that this “open-loop” method is only an approximation. The flexible 
fuselage dynamics are not part of the coupled rotor-fuselage free flight trim procedure. The 
calculated hub load does not include the flexible fuselage motions. 


a y 

a z 

ai 




-I cockpit, 4/rev 


N x N y 




Fy 

F z 

M x 

My 

M z 


- hub, 4/rev 


(2.20) 


Figures 2.4a and 2.4b show the 4/rev acceleration magnitudes at the pilot and copilot 
stations. Note that the result from HeliUM closely follows the UMARC result up to 90 
kts. Beyond 90 kts, the HeliUM result continues to rise as the forward speed increases. 
Overall the predictions qualitatively follow the trends of the flight test data except in a 
higher speed range. This over-prediction could be caused by several factors. First, the 
hub load calculation from HeliUM does not include the effect of the flexible fuselage 
dynamics. The effect of aerodynamic damping on the hub load calculation is also not 
considered. Second, HeliUM does not have any passive vibration damping device such as 
the hub absorbers, the bifilars, and the spring-mass fuselage absorbers. Third, equation 2.20 
assumes that the 4/rev cockpit station acceleration is a linear combination of the 4/rev hub 
shears and moments. Because vibration is not a linear phenomenon, this assumption may 
not hold true in the high-speed flight condition. The cockpit station accelerations provided 
here are only intended to serve as a qualitative measure. 


2.8 Optimization formulation 

In the first attempt to formulate the optimization problem, the trim equations were included 
directly in the form of equality constraints // (X) (recall that the trim problem is formulated 
as a set of nonlinear algebraic equations as stated in section 2.4.1, and the trim unknowns 
X were included as design variables. Therefore, the optimization problem was formulated 
as follows: minimize the norm of 4/rev in-plane hub shears, F\ p 

F(X) = H-F 4 PH 2 — > min 

Subject to 

Equality Constraints, hj(X) < e 

Of the 29 equality constraints, 11 represented trim conditions for the entire aircraft, 4 for 
the inflow trim equations, and 14 for the main rotor equations. The vector X of design 
variables was composed of 31 elements, namely, 29 trim variables, and the sine and cosine 
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magnitudes of the HHC input. The initial solution was obtained from an algebraic trim 
procedure without a HHC input, and therefore it was always feasible. The optimization was 
carried out using a modified method of feasible directions (ref. 55), as implemented in the 
code DOT (ref. 56). The numerical properties of this formulation proved to be extremely 
poor. Convergence was very slow, and the software often terminated the optimization for 
lack of progress. Several variations of the baseline process were tried unsuccessfully and 
this formulation was abandoned. 

A different approach to the optimization process proved more successful. The problem 
was formulated as an unconstrained minimization: 

F(X) = ||F 4 p|| 2 — » ► min 

Subject to 

Unconstrained optimization 

with a vector X of design variables consisting of just 2 elements, namely the sine and 
cosine magnitudes of the HHC input. This way, the trim procedure is decoupled from the 
optimization, and it is simply executed separately for every value of X proposed by the 
optimizer. The optimization was carried out using a Broyden-Fletcher-Goldfarb-Shanno 
(BFGS) algorithm (ref. 55), as implemented in the code DOT. 
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Figure 2.1. Cockpit vibration comparison; 18,000 lb, linear inflow model, rigid fuselage. 
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Figure 2.2. Cockpit vibration comparison; 18,000 lb, free wake model, rigid fuselage. 
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Figure 2.3. SH-60 fuselage NASTRAN model; courtesy of M. Yang and I. Chopra. 
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Figure 2.4. Cockpit vibration comparison; 18,000 lb, linear model, flexible fuselage. 
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3 Active Rotor Control System for Vibration Suppression 

The active rotor control system is implemented in the nonlinear helicopter simulation as 
shown in figure 3.1. This loop consists of three main parts: the harmonic analyzer, the HHC 
controller, and the discrete HHC update. Because the final higher harmonic control (HHC) 
and automatic flight control systems (AFCS) interaction study (chapter 5) is performed 
in the continuous linear time-invariant system, each component in the feedback path is 
converted to an equivalent linear model. 

This chapter is divided into three main parts. The first describes the function of the 
harmonic analyzer and its linear time-invariant equivalent model. The second describes 
the algorithm of an HHC controller for the steady-state vibration suppression. The third 
section describes the discrete HHC update and its linear time-invariant equivalent system. 

3.1 Harmonic analyzer 

In several research studies, the method of extracting n/rev vibration components is to use 
a harmonic analyzer. The harmonic analyzer can be formulated by using either an analog 
bandpass filter (refs. 6, 9, 32) or a Fourier analyzer (refs. 7, 10, 29, 30, 57). 

3.1.1 Analog bandpass filter method 

This type of harmonic analyzer consists of three components: a bandpass filter, 

demodulator, and a lowpass filter (fig. 3.2a). An analog filter operates on continuous-time 
signals, and provides a continuous sensor output without the effect of the sampling window 
(Sec. 3.1.3) that is typically associated with Fourier analysis. 

First, the analog bandpass filter extracts the spectral band of interest from the source 
signal. This spectral band is centered on 4/rev frequency and has a width of c u B w- The 
pre-filtered signal is demodulated by multiplying the exact 4/rev harmonic frequencies. 
The resultant signal contains all the sum and difference frequencies created by the 
multiplication. Finally, the lowpass filter removes frequency above ujbw/2 with an 
assumption to < lu bw /2. The following example illustrates this process. Let Z 4P (t) be 
one of the spectrum band of some general non-periodic hub load signal Z(t). 

Z 4P (t) = At cos(4f2t + cut) (3.1) 

where A 4 is the 4/rev amplitude, H is the rotor speed in radian per second, and tut is the 
4/rev phase angle in radian. Although the hub loads in a trim condition are restricted to 
periodic waveforms, there is no such restriction during gusts and maneuvers in which the 
hub loads may contain significant non-periodic transients. Therefore, the source signal 
Z (t) is first screened through the bandpass filter to extract its components near the 4/rev 
frequency, i.e., more precisely, in the range of 4Q — tu BW /2 an d 4Q + lu bw /2. This pre- 
filtered signal is called Z 4P (t) and can be written in the form of equation 3.1. Next, the 
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pre-filtered signal is demodulated, i.e., multiplied by cos4fl/ and sin 40/ as follows: 

A 4c = Z 4P (t) cos 40/ 

= -A 4 [cos (cut) + cos(80/ + cut)] (3.2) 

A 4s = Z 4 p(t) sin 40/ 

= -A 4 [— sin(cut) + sin(80t + cut)] (3.3) 

The demodulated signals A 4c and A 4s are passed through the lowpass filter to remove all 
frequencies above lubw/2 and doubled with an assumption cu < ubw/2- The resultant 
signal is given by equations 3.4 and 3.5. 

A 4c = A 4 cos(cut) 

A 4s = —A 4 sin (cut) 

Equation 3.1 can be rewritten, using equations 3.4 and 3.5, as 

Z 4 p (t) = A 4 cos (40/ + cut) 

= A 4 cos(cut) cos(40/) — A 4 sin(cut) sin(40/) 

= A 4c cos(40/) + A 4s sin(40/) (3.6) 

Using analog bandpass filter to extract the 4/rev signal adds a large time delay to the 
system because the method requires a high order bandpass filter with narrow passband 
width (small c opw) 10 extract the steady-state vibration value. The analog harmonic 
analyzer does not present a problem for steady-state vibration extraction, however the large 
time delay will mask all the transient responses. 

3.1.2 Fourier analyzer method 

Another method of extracting the harmonic components from the source signal is to use a 
Fourier analyzer, which was applied in this study. This type of harmonic analyzer consists 
of three components: a sample window, the Fourier analyzer, and a lowpass filter (fig. 3.2b). 
The sample window serves as the data buffer which stores incoming data streams. The 
Fourier analyzer then identifies the harmonic contents of the source signal within the 
sample window. This Fourier analyzer can be either a Fourier series approximation or a 
Fourier transform in either the continuous-time or discrete-time domain. The lowpass filter 
then removes the undesired frequency contents above n/rev frequency. The Fourier series 
approximation method was chosen as the Fourier analyzer within the harmonic analyzer 
because the HHC/AFCS interaction study is performed in continuous-time domain system. 
Nevertheless, one can use the Fourier transform (such as FFT, Appendix B) as the Fourier 
analyzer for the digital version of the harmonic analyzer. 


(3.4) 

(3.5) 
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The theory of Fourier series approximation lies in the idea that most signals, and all 
engineering signals, can be represented as a sum of sine waves: 

^ OO 

/W = + X! [a n cos(27rn/f) + b n sm(2nnft)] 

^ n = 1 


with 


a n = ^-[ f(t)cos(2nnft)dt, n = 0,1,2,... 
1 Jo 

b n = ^ [ f{t) sin(27 rnft) clt, n = 1,2,3,... 
1 Jo 


(3.7) 


where T is the fundamental period and / = 1/T is the fundamental frequency in Hz. For 
example, the vertical vibratory hub load Fz can be approximated using the finite version of 


equation 3.7 that will pass through N data values in one fundamental period: 

x / 2 /2-7T nk\ ( N/2)-i f 2Trnk\ 

F z (kAt) = Fy+£F Znc cos(^J+ £ F ^(^f) (3.8) 

with 

t = kAt, k = 1,2 ,N, where At = T /N (3.9) 

1 n 

Fz = ^J2 F z(kAt) (3.10) 

k = 1 

9 N ( r >'Trnk\ N 

F Znc= J2F z (kAt)cos — j , n= 1,2,..., — -1 (3.11) 

F Zns = J2 F z( kAt ) sin (“77“) > n= l>2,...,y-l (3.12) 


The fundamental frequency / is the rotor speed f! in rad/sec or Q/2 tt in Hz, and the 
fundamental period T is 27r/f2 second. To extract nl rev components of F z , the sampling 
frequency must be at least twice as fast as nl rev frequency to avoid aliasing problems. In 
this study, a factor of 6 is chosen which leads to a sampling frequency of (6r2.fi/27r) Hz. 

For a rotor with four identical blades and zero tracking error 1 , the only frequencies 
transmitted to the fixed system are the four multiples per revolution (4/rev, 8/rev, 
12/rev . . .). Therefore, the sampling frequency required to extract F Z4:C and F Z4s of the 
Sikorsky UH-60 helicopter with a nominal rotor speed fi=27 rad/sec is f s = 6nfi /2 tt = 
(6 x 4 x 27)/27 t = 103 Hz or N = 24 sample data per rotor revolution. Although the 


'The main rotor blades are all flying in the same tip-path-plane and maintain equidistant angular spacings 
during flight. 
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major vibratory hub loads interested in this study have a 4/rev frequency, the study also 
monitored 8 and 12/rev frequency, which are the second and third harmonics transmitted to 
the fixed system for a four-bladed helicopter. In this case, the sampling frequency required 
is f s = 6nQ/2n = (6 x 12 x 27)/27t = 310 Hz or IV = 72 sample data per rotor revolution. 

The calculation of the Fourier coefficients is very time consuming since it requires N 2 
number of function evaluations. To reduce the computation time, the most frequently used 
algorithm for real-time applications is the fast Fourier transform (FFT, Appendix B). It is a 
discrete Fourier transform that reduces the number function evaluation from N 2 to N log N. 
Since the helicopter simulation program used in this study is in the continuous-time domain 
and it is not a real-time simulation, the simulation time does not advance to the time frame 
until the Fourier series calculations is finished. Therefore, the additional computation time 
required for the Fourier series calculations has no impact on 4/rev vibration extraction. 

The lowpass filter implemented in this harmonic analyzer is a 4 th order Bessel lowpass 
filter with the break frequency tu 0 at 6.5/rev. The 6.5/rev break frequency is chosen to 
produce a -12 dB magnitude drop between 4/rev and 8/rev signals. Additional information 
about Bessel filter will be discussed in section 3.1.4. 

Use of the Fourier analyzer (either Fourier series approximation or FFT method) to 
extract the 4/rev signals causes additional delay to the system. The source of delay is 
from the sample window, which is discussed next. 


3.1.3 Effect of windowing 

When performing a digital harmonic analysis with a physical system, a sample window 
must be used, as it is necessary to truncate long data streams to a finite size. The size of the 
window has a significant effect on the accuracy of the extraction of the desired frequency 
components. A large window; i.e., a window that extends over a long time, increases the 
accuracy of the low-frequency components identification but degrades the high-frequency 
identification. On the other hand, a small window improves high-frequency components 
identification but degrades the low-frequency. Generally, the minimum window size is one 
cycle of the source signal. For the 4/rev hub load study, the minimum sample window is 
equal to quarter revolution. However, a sample window of one revolution (4 cycles of the 
source signal) was used in this study. As the rotor rotates beyond its first revolution, the 
sample window advances with it continuously. 

Figure 3.3 illustrates the time delay introduced by the sample window. The vertical 
hub load Fz, shown in the second figure, starts from a trim condition without HHC input 
for the first two revolutions. At the end of the second revolution, a 4/rev HHC input is 
added (this is an arbitrary input, which will not necessarily reduce vibrations), and the 
helicopter reaches the new steady state condition. In the third revolution, Fz has reached 
the steady state almost instantaneously. Although there is a low frequency drift, mainly 
1/rev response, the third revolution is dominated by the 4/rev response, but is near to the 
new steady state condition. The spectral analysis performed on the third revolution also 
confirms this finding and the result is shown in the third figure. However, according to 
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on-line Fourier analysis with a moving sample window (fourth figure), F z 4C and F z iS take 
approximately one rotor revolution to reach the new steady state condition. This does not 
agree with the result stated earlier. The cause of this difference is the sample window. In 
other words, the sample window behaves as a lowpass filter, which adds large time delay 
and masks all transient responses. 

3.1.4 Equivalent lowpass filter 

A window essentially behaves as a lowpass filter. The sample window used in the study is 
based on a rectangular or ’’box car” window. Figure 3.4a is the rectangular window, h(t), 
in time domain which has a window size of 2T 0 . Its expression is given by: 


where / is the frequency in Hz. Figure 3.4b shows that the Fourier transformation of a 
rectangular waveform consists of a central lobe which contains most of the energy of the 
window and the side lobes which generally decay rapidly. The magnitude difference of the 
first two lobes is 13.4 dB (79% reduction) with a break frequency of 1/2T 0 Hz. 

Equation 3.14 is a closed form solution, and is a function of frequency. The LTI system 
analysis, equation 3.14 can be approximated by an equivalent lowpass filter. The equivalent 
lowpass filter chosen is the Bessel filter because it has the following characteristics: 

• k poles and no zeros 

• DC gain = 1 

• Break frequency = u 0 

• Maximally flat group delay about 0 Hz, and the phase response is approximately 
linear in the passband 

• The linearity degrades at the higher frequencies, and the group delay drops to zero 

• No overshoot around break frequency 



1*1 < T 0 
|*| > T 0 


(3.13) 


and its Fourier transform is given by: 



(3.14) 
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For the Sikorsky UH-60 helicopter, the Bessel filter chosen is a 4 th order function with the 
break frequency lu c at 1/2 T 0 . 2 T 0 is the length of the sample window, which is equal to the 
time to complete one rotor revolution. 

The frequency response comparison II (/) of the rectangular window and that of the 
Bessel filter are compared in figure 3.4 (c). The Bessel filter low-passes the signal at a 
break frequency oj 0 of 4.3 Hz and produces a frequency drop off similar to the rectangular 
sample window. Note that the Bessel filter designed in this section is only implemented in 
an LTI system analysis to mimic the dynamics and the time delay of the actual rectangular 
sample window. 


3.2 Higher harmonic control algorithm 

3.2.1 T-matrix method 


The closed-loop HHC algorithm implemented is based on the fixed-gain T -matrix feedback 
controller: 

Z 4 p (k) = Z 4P (k — l) + T [G hhc (k) - G hhc (k - 1)] (3.15) 

Equation 3.15 is a difference equation for discrete-time domain system. The variable k is 
the discrete-time index, while Z 4 p is the vibration response vector consisting of cosine and 
sine components of 4/rev vibratory hub loads excluding the 4/rev yawing moments: 

Z 4 p = \ F x 4C i Fx 4 g i Fy 4C , Fy 4S , Fz iC i Fz 4S , Mx 4C i Mx 4S , My 4C , My 4S ] T (3.16) 

and is a function of the state vector x, the pilot inputs G p u ot , and the HHC inputs Ghhc 


Z 4 p — 

f(x, 0 pilot 1 ) @hhc) 

(3.17) 

L ) 

pilot 

$lom 3 coh 

(3.18) 

@hhc 

[# 3 C, # 3 Si G&C i Oa S, O5C, 65 s] T 

(3.19) 


The T-matrix is the Jacobian of function f computed about a reference input vector, G hhCo '- 


T 


dG 


G 


hhco 


(3.20) 


In other words, T is a linear approximation of the 4/rev vibration response Z 4 p to the HHC 
inputs G hhc at a steady-state condition. That is, equation 3.20 assumes that changes in the 
vibration response AZ 4P with respect to the changes in the HHC input A G hhc are linear 
over the entire range of Ghhc ■ This relationship can be written as 


AZ 4 p — T A G hhc 


(3.21) 


In this study, the helicopter is trimmed without the HHC input; therefore, the reference 
input vector G hhCo is a zero vector, and the A G hhc in equation 3.21 is the same as G hhc . 
Total 4/rev vibration approximated using T -matrix method is given by 

Z 4P = Z 4 p 0 + T Ghhc (3.22) 
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where Z 4 p 0 is the 4/rev vibrations of the nonlinear baseline (HHC-off) case. 

For vibration suppression, optimal control is obtained by minimizing the cost function 
J: 

J=\ Z lp(k) Q z 4 p(fc) + l - e T hhc {k) R e hhc (k) (3.23) 

where 0 and R are the weighting matrices on the responses and controls: 


Q = diag{l, 1, 1, 1, 1, 1, q 7 , . . . , gio} (3.24) 

T = diac/{l, 1, 1, 1, 1, 1} (3.25) 


and q 7 , ... ,qio = 1/A z 2 cg where A z cg is the vertical displacement of the rotor hub to the 
center of gravity of the helicopter. The choice of the weighting 1/ /S.z 2 cg transforms the 
moments to the equivalent forces. The optimal control is computed by setting the first 
derivative of the cost function of equation 3.23 to zero and solving for the optimal HHC 
input: 



(3.26) 


With this scheme, the HHC input is computed based on the current response vector: 


9{k) = T f T t 0{k - 1) - T f Z 4P (fc - 1) (3.27) 


where the fixed-gain regulator is 

T f = (T t QT + R)' 1 T t Q (3.28) 

If R = 0 or T t QT R, T f is a pseudo-inverse of T, and equation 3.27 becomes 

0(k) = 0{k - 1) - T f Z 4 p(/c - 1) (3.29) 


3.2.2 T-matrix validation 

As stated before, T -matrix is a linear approximation of the 4/rev vibration response Z 4 p to 
the HHC inputs 0 hhc at a steady-state condition. The total 4/rev vibration response Z 4 p of 
the nonlinear model and that of the T-matrix approximation are compared in figures 3.5- 
3.7 for 3, 4, and 5/rev inputs to determine the accuracy of the T-matrix approximation. 
The total 4/rev vibration response Z 4 p of the nonlinear model at the steady-state condition 
was computed for HHC input, and the 4/rev vibration response was extracted from the 
helicopter hub loads using Fourier series approximation. The total 4/rev vibration response 
Z 4 p of the T-matrix approximation is computed using equation 3.22 where Z 4 p 0 is the 
4/rev vibrations of the baseline (HHC-off) case from the nonlinear simulation. The HHC 
input for both methods have an amplitude of 1° with a phase angle varying from 0° to 360° 
with increment of 30°. 

The diamond symbol represents the baseline (HHC-off) 4/rev vibration responses Z 4 p 0 
from the nonlinear model, with values tabulated in table 3.1. The open circles represent the 
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vibration responses from the nonlinear model with the HHC inputs engaged; solid circles 
represent the vibration responses from T-matrix approximation based on equation 3.22 
with 6 hhc determined from equation 3.29. The number next to the symbol is the n/rev input 
phase angle. 

These figures illustrate the 4/rev vibration prediction error resulting from the T-matrix 
approximation. In the 3/rev case (fig. 3.5), 4/rev vibration responses from the T-matrix 
approximation match well with that from the nonlinear model. For the 4/rev and 5/rev 
cases, there are differences between the two methods. Those difference are from an earlier 
assumption that the vibration response to the HHC input is linear over the entire range of 
Ohhc- Since the vibration responses to the HHC inputs are not necessarily linear, the small 
differences between linear and nonlinear models are expected. 

3.3 Discrete HHC update 

The ideal HHC inputs computed by the T -matrix controller for vibration suppression are 
not returned to the rotor system at every time step. The HHC input has a discrete update rate 
which typically varies from 0.5 to 16 times per rotor revolution (refs. 5-10). A typical HHC 
input update rate is once-per-revolution. In the discrete-time domain, the discrete HHC 
update is performed by the sample-and-hold operation. To implement this in a continuous- 
time domain system, the effect of the sample-and-hold operation must be approximated in 
the continuous-time domain. 

Figure 3.8 illustrates the effect of a sample-and-hold operation on a continuous signal. 
The sampler transforms the continuous signal to an amplitude-modulated pulse signal at a 
sample rate u s . At the output of the digital controller, the digital signal must be converted 
to analog by the process called digital-to-analog conversion. The simplest device that 
transforms digital input to analog output is a zero-order-hold. The bottom of figure 3.8 
shows the relationship between digital input and analog output. The zero-order-hold holds 
the value of the sampled signal over T s second to produce a reconstructed signal with 
staircase waveform. Notice that an approximation to the reconstructed signal is identical 
to the original signal with a delay of T s /2 second. Therefore, a zero-order-hold operating 
at a sample rate tu s is equivalent to a time delay of tt/lu s second. Similarly, the discrete 
HHC input operating at cu s2 frequency also can be approximated by a Pade function with a 
time delay of n/cu S 2 second. 
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Table 3.1. Baseline (HHC-off) vibration level. 



4P Cos-Comp. 

4P Sin-Comp. 

Amplitude 

F x (lb) 

151.6 

87.8 

175.2 

Fy (lb) 

73.5 

-61.3 

95.7 

F z (lb) 

39.5 

8.9 

40.5 

M x (ft-lb) 

40.1 

62.6 

74.3 

My (ft-lb) 

80.0 

30.2 

85.5 
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Figure 3.1. Nonlinear simulation scheme; multi-rate system. 
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Figure 3.2. Internal components of the harmonic analyzer. 
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Figure 3.3. Example of the time delay induced by the sample window. 
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Figure 3.4. Frequency response of the rectangular window. 
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Figure 3.5. 4/rev vibration response comparison; 2L 3 =1°; V=120 kts (/x= 0.28), 
Weight= 14,000 lb. 
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Figure 3.6. 4/rev vibration response comparison; At=l°; V=120 kts (/i=0.28), 
Weight= 14,000 lb. 
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♦ Baseline case 
O Nonlinear model 
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Figure 3.7. 4/rev vibration response comparison; 2L 5 =1°; V=120 kts (/x= 0.28), 
Weight= 14,000 lb. 
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Continuous Signal 





Figure 3.8. Time delay approximation for sample and hold circuit (sampler and zero-order- 
hold operating at a same rate, c u s . 
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4 Extraction of the Constant- Coefficient Linearized 

Model 


A key ingredient for the study of potential interactions between HHC and flight control 
system is a linearized time-invariant model of the helicopter dynamics, including higher 
harmonic inputs and controls. This chapter contains the derivation of such a model, and is 
composed of three sections. 

Section 1 summarizes the main steps of the extraction of a conventional linearized model, 
i.e., one without higher harmonic inputs and controls. Section 2 extends the derivation 
to include such higher harmonics to show that: (i) one portion of the output equation 
is the equivalent of the traditional T-matrix, and (ii) through an appropriate formulation 
of the output equation, the need for online identification and adaptation of the T-matrix 
in maneuvering flight is substantially reduced. Section 3 describes the application of 
the methodology to simple linear rotor equations, for which analytic expressions for the 
coefficients of the model can be derived. 


4.1 Extraction of a linearized model without higher harmonics 


Consider the equations of motion of the helicopter written in symbolic form as: 

f(x,x,u;V>) = 0 

and take first order differentials 

df(x, x, u; ijj) — 0 

which can be expanded into 

Of 

Ox 


x = x 0 


Of 

dx + Tvr 
ox 


, Of 

dx + 7 ^ 

OU 


du = 0 


(4.1) 

(4.2) 

(4.3) 


u = u 0 


x = x 0 

where the subscript (. . ,) 0 denotes the trim values of the respective vectors. Replace now 
d(. . .) with A(. . .) and introduce the notation 


[£«0] 

def 

Of 

Ox 


(4.4) 



X = x 0 


Hi M] 

def 

Of 

Ox 

X = X 0 

(4.5) 

[BiWO] 

def 

Of 

Ou 

U = Uq 

(4.6) 


Then equation 4.3 can be rewritten as: 

Ak = - [-^(VOr 1 [A (VO] Ax - [E{^)}~ 1 [B 1 (i/j)\ Au 

= [d(^)]Ax + [B(i))]Au (4.7) 
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with [A(ip)\ d = -[E(ip)\ and [B(ip)] '= ~[E(ip)] 1 [B^ip)]. The 

linearized matrices [E(ip)], [A^ (ip)], and \B\ (ip)] can be calculate using finite difference 
approximations. For example, using central finite differences, the j-th columns of the 
matrices [Ai(ip)], [B i(ip)], and [ E( ip)] at the azimuth^ are given by, respectively: 




{Bityi)}, 

{E(m j 


di 

dxj 

df 

duj 

<9f 

din 


X = X 0 


U = U 0 


X = x 0 


f (x 0 + hej] 1/jj ) - f (x 0 - he j ]'ijji) 
2 h 

f (u 0 + he^ipj) - f (u 0 - hej; ipj) 
2 h 

f (x 0 + hej] 'i/jj ) - f (x 0 - 

2 h 


(4.8) 

(4.9) 
(4.10) 


where e ? is a vector with all its elements equal to zero except for the j-th, which is equal 

to one, and h is the finite difference step size. All the matrices above are periodic, with 

common period equal to one rotor revolution. Therefore, the state matrix [A(r/>)] and the 
control matrix [B(iji)\ are also periodic, and can be expanded in Fourier Series: 

K 

[A(ip)\ = [A 0 ] + ( [ A kc] cos + [^fes] sin k^) (4. 1 1) 

k=\ 

K 

[B(ip)} = [5 0 ] + J2 (i B kc\ cos kip + [B ks ] sin kip) (4.12) 

k = 1 


If the state vector x is defined entirely in a fixed coordinate system, then a time invariant 
linearized model can be obtained by retaining only the constant matrices [A 0 ] and | /i 0 1 . If, 
additionally, the blades are assumed to be identical, then the summations in equations 4.11 
and 4.12 only contain harmonics that are multiples of the number of blades. Therefore, for 
an A r -bladed rotor, k = N,2N,3N, . . . 


4.2 Extraction of a linearized model with higher harmonics 

This section presents the extension of the linearization procedure to the case in which 
both the state vector x and the control vector u contain higher harmonics. The precise 
definitions of x and u will be introduced first, together with general expressions for the 
linearized system. The derivation of the control matrix [B{ip)\ will be presented next, 
as it requires only minor modifications of the baseline procedure of the previous section. 
Finally, the derivation of the state matrix [A{ip)\, which requires some special treatment, 
will be presented. 
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4.2.1 Definitions 


The control vector u(r/>) used in the present study is defined as: 

u(V0 = 

where u p ,i 0 i (vj) is the vector of conventional pilot controls 


Hpilot (VO 

U HHCW 


(4.13) 


Hpiloti^') ^ lot, 3 Ion $col 

and u HH c is the vector of higher harmonic controls 


Spedf 

(4.14) 

Obc $5s] T 

(4.15) 


The HHC is applied to the blade in the rotating system. In the 4-bladed helicopter 
configurations used in this research, 3/, 4/, and 5/rev control inputs in the rotating system 
are required to generate the desired 4/rev inputs in the fixed system. The vector u(/0) should 
be interpreted as “perturbations from the trim values of the controls”. The state vector x(C), 
also representing perturbations from trim values, can be written in the symbolic form: 


x(V0 


X£ 

X MR 


where xp is the vector of states not associated with the main rotor, defined as: 


xp = [u v w p q r </> 6 -0 A 0 A c A s \ tr u x u y ] T 


(4.16) 


(4.17) 


and ~x-mr is the vector of rotor states, defined in a fixed coordinate system. The elements of 
the rotor state vector are based on the assumption that each state is composed of an average 
and a 4/rev portion, both azimuth dependent. For example, with the longitudinal rigid body 
flapping /3i c (-0) written as: 

PiM = PicM) + /?1C4C WO cos4 ^ + Pi c 4 s (V 0 sin Aip (4.18) 

the quantities /?i Ca „ e ('0), Pi C4c ('0)- an d fi\ Cis (VO will be considered as states and included 
in the rotor portion x A/ r of the state vector. The state T lc „, ue ('C) is equivalent to the 
longitudinal flap state that would appear in a traditional rotor state vector. The additional 
higher harmonic states T 1C4r: (C) and {%)) represent a new way of modeling the effects 
of higher harmonic control, introduced for the first time in the present research. Although 
the formulation of equation 4.18 appears intuitively reasonable, it will not be justified on a 
rigorous theoretical basis. However, its validity will be established through simulation, by 
comparing linearized and nonlinear responses to pilot inputs. 

The assumption that each rotor state is composed of an average and a 4/rev portion leads 
to an expanded state vector defined as follows: 


x = 


X 4P 


(4.19) 
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where x,„ ;e contains the vector x B defined in equation 4.17 and the average rotor states, 
that is: 


x 5 0o ave C-ave Save fi^ave 0o ave /^1 Cave Save /^2 a 

* * * f^OaTje ^lCa^e (*4s a „ e C°2aue COaue Cf Save C2ave 

PlCave 01 S 


ICave t is. 


aue / ^ave / 'Jave r i^ave r ±a ave / ^ave 


and x 4 p contains the 4/rev components, sine and cosine, of the rotor states: 


x 4P = 


0 O 4 C 0 O 4 S 01 C 4 c 01 C 4 S 01 S 4 c 01 S 4 s 024 c 024 s ■ ■ 

■ ■ ■ 00 4c 0o 4s /^lC4c /^lC4s /^lS4c /~^lS4s /^24 c /~^ 24 s 

* * * C^04c Co 4 s C"lC 4 c C^lC 4 s C^lS 4 c ClS 4 s C^ 4 c ^24 s • • 

Co 4 c Co 4 s ClC 4 c ClC 4 s olS 4 c jlS 4 s j 24 c j 24 s 

■ ■ ■ 0O 4c 0O 4 s 01C4c 01c 4s 01«4c 01«4s 02 4 c 02 4a 

• • • 004c 004s 01C4 C 01C4s 01 *4c 01s 4 s 02 4c 02 4s 


(4.20) 


(4.21) 


The notation in equations 4.20 and 4.21 reflects the fact that in the present study the rotor 
blades are modeled using one rigid flap mode /3, one rigid lag mode and one flexible 
torsion mode </>, but both equations can be rewritten for a generic number of rigid and 
flexible modes. Also note that both vectors x, n , e and x.hhc are i n general time dependent. 

With these definitions of the state and the control vector the linearized system [Eq. 4.7] 
becomes 


I -X - ave 

L=[ 

{ * 4P 

[ [ 


r *-ave 

-4-21 


-4-12 

Ahhc 


X-ave 

X 4P 


Bave B\2 

B 21 b hhc 


U pilot 
U HHC 


(4.22) 


where now all the partitions of A and B are time-invariant. In other words, by decomposing 
the state vector into an average and a 4/rev component, the original linearized system with 
periodic coefficients has been converted into a larger linearized system, but with constant 
coefficients. 

The linearized model also includes an output equation, which has the form: 


f y ave 


' / 

0 ' 

1 y 4 P 


0 

C 22 

I F ave 

> — 

cv, 

C 32 

{ F 4P J 


. c 41 

C 42 . 



def 

= y 


dej q 


0 0 

0 0 

-D 31 -D 32 

-D 41 D42 

def jj 


H pilot 

UHHC 


(4.23) 


where the C and D matrices have constant coefficients. The vectors F„,, e and F 4 p contain 
average and 4/rev hub loads at the hub, and are defined as: 


F a ve — [F Xave F Xave Fj 
F ip = [F 


V M x 
yave ^ave 


My ave M Zave 


X4c F X 4s Fy 4c Fy 4s F Zic F Z4s M X4c M Xis My 4c My 4s M Z4c M Z4s 


(4.24) 


(4.25) 
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where F x , F y . F z , and M x , M y . M z denote the rotor force and moment components along 
and about the body axes. The remaining two partitions of the output vector y in 
equation 4.23 are y ave and y 4P . The output subvector y ave is identical to the average state 
vector x ave ■ The output subsector y 4P is the global 4/rev rotor state vector: 


Yap 


A) 4c A) 4s 01 C4c 01 C 4s 01s Ac 

c' c' c' c' c' 

• * * S04c ^04s S lc4c j lC4s SIS4 


0'lSAs 02 4c 02 As ' ' ' 

C C C 

c V>5 lS4s j24c ^24s 



(4.26) 


The portion of the output equation corresponding to y ave and y 4P is simply a mathematical 
means to indicate that the outputs are the average and global 4/rev rotor states; no physics 
are involved. 

The submatrices C 31 and C 32 express a linearized relationship of the average hub loads 
F ave with the average rotor state x ave and the 4/rev rotor states x 4P . Similarly, the 
submatrices C 4 i and C 42 express a linearized relationship of the vibratory loads F 4P with 
the average rotor states and the 4/rev rotor states x 4P . 

As for the feedforward matrix, the submatrices D 31 and D 32 link the average vibratory 
loads to pilot and HHC inputs, respectively. The submatrices D 41 and D 42 link the 4/rev 
harmonics of the vibratory loads to pilot and HHC inputs, respectively. Therefore, the D 42 
submatrix is the equivalent of the T-matrix in typical HHC studies. The submatrix D 41 
represents the effects of pilot maneuvers on the vibratory loads: these effects are not taken 
into account explicitly in typical HHC studies, instead, the maneuver effects are captured 
by online identification of the T-matrix and adaptation. By including the maneuver effects 
in the output model, the need for adaptation is substantially reduced. 


4.2.2 Extraction of the control matrix B 

The extraction of the control matrix B is presented first, because the procedure is more 
similar to that for the traditional linearization without higher harmonic components of the 
states. In fact, the control perturbation vector u is already defined in the rotating system. 
The control matrix B is extracted through numerical perturbation of the full nonlinear 
equations of motion about a trimmed equilibrium position. Each element of the matrix 
is obtained using central difference approximations. The calculation proceeds as follows: 
For every azimuth angle 


1. Perturb the k - th element Uk('iPi) of the control vector u (pilot and HHC controls are 
treated in exactly the same way) by A u k , i.e., let the perturbed control vector u + (T l ) 
be 
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U 1 
U2 



(4.27) 


where the subscript “+” denotes the positive perturbation in the central difference 
calculation. 

2. Substitute the perturbed control vector u + (00 into the system of equations of motion 
of the helicopter, to obtain the perturbed acceleration vector x/,> + 


where a subscript R has been added to the state vector to indicate that the rotor 
portions are formulated in the rotating system (note that the state vector in the 
linearized model is entirely expressed in the fixed system). The state vector 
x/,. corresponds to the desired trim condition, and is held constant during the 
perturbation. 

3. Repeat the two previous steps with a negative perturbation of the k - th control, Uk — 
Auk, to obtain the perturbed acceleration vector xr_. 

4. Build the derivative using central difference approximations. This derivative is the 
k - th column of the B R matrix (i.e., with the rotor portions still in the rotating system) 
at the azimuth angle 0*: 


5. Repeat the four previous steps for each of the m elements of the control vector, i.e., 
for Uk, k = 1, . . . , m, to obtain the complete control matrix 


The next step of the linearization procedure typically consists of performing a multiblade 
coordinate transformation, to convert the rotor states from the rotating to the fixed system, 
and therefore to obtain a control matrix 5(0*) entirely in the fixed system. After steps 
1-5 are made for a sufficient number of azimuth angles 0$, the resulting control matrices 
5(00 are typically averaged to obtain the final constant control matrix B. This is the 
traditional linearization procedure used in the present study for the calculation of the rows 
of the B matrix corresponding to the “average” states, i.e., for the submatrices B ave and 
B 12 in equation 4.22. 


± R+ = f(x fl ,U +,0i) 


(4.28) 



(4.29) 


BrM = {Br(A)} 2 ■ • • {Br(A)}J 


(4.30) 
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Some additional manipulations are required for the rows corresponding to the 4/rev 
states, i.e., for the submatrices B 2 1 and B HHC . These manipulations consist of Fourier 
analysis of B(ipi) to extract the 4/rev cosine and sine harmonics. Define: 




B 4c = — B (i/>i) cos 4 

iV V> i=i 

(4.31) 



B 4s = -T sin 

^ i= 1 

(4.32) 


where is number of azimuth angle in one rotor revolution. Then it is essentially 


[£>21 b hhc \ 


B 4c 

B^ 


(4.33) 


except that the rows of B 4r and B 4s must be appropriately permutated because the state 
sub vector x 4 p [Eq. 4.21] is arranged with the 4/rev cosine and sine components interlaced 
rather than grouped together. 


4.2.3 Extraction of the state matrix A 

The general procedure to extract the state matrix A is similar to that of the control matrix 
B, except that the state vector is defined in the fixed system, both for the average and the 
4/rev components. 


4.2.3.1 Rows corresponding to the average states x,„, e 


The rows of the A matrix corresponding to the average states x„ w , , i.e., the submatrices A ave 
and A 2 i in equation 4.22 can be obtained with the same procedure as previously shown for 
the B matrix, i.e., through the following steps. 

For every azimuth angle Vv 


1 . Perturb the k - th element x avek ('fy) of the partition x„,, e of the state vector x by Axk, 
i.e., let the perturbed state vector x + (?/y) be 


%ave i 
%ave 2 


X+O*) = < 


•Eave^ > 


%ave]y 


X 4 p 


(4.34) 


where the subscript “+” denotes the positive perturbation in the central difference 
calculation. 
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2. Substitute the perturbed state vector x + (00 into the system of equations of motion 
of the helicopter to obtain the perturbed acceleration vector x + . Because the rotor 
equations are formulated and implemented in the rotating system, x + (00 must first 
be converted to the rotating system, using a multiblade coordinate transformation 
that yields the corresponding rotating state vector x R+ (00 


where the subscript R again indicates that the rotor portions are in the rotating 
system. The control vector u corresponds to the desired trim condition, and is held 
constant during the perturbation. 

3. Repeat the two previous steps with a negative perturbation of the k - th average state, 
x aV e k ~ A-Xk, to obtain the perturbed acceleration vector x ft _. 

4. Build the derivative using central difference approximations. This derivative is the 
k - th column of a matrix Ap(00 at the azimuth angle 0; , that is: 


Position in the state matrix and dimensions of 2Lp(00 are the same as the submatrix 
A ave in equation 4.22, but A ave is constant and in the fixed system, whereas A R (00 
is periodic and in the rotating system. 

5. Repeat the four previous steps for each of the N elements of the state vector partition 
x a „ e , i.e., for x avek , k = 1, . . . , N, to obtain the complete matrix Ar(00 


The next step of the linearization procedure typically consists of performing a multiblade 
coordinate transformation, to convert the rotor states from the rotating to the fixed system, 
and therefore to obtain a state matrix A (00 entirely in the fixed system. Then, after 
steps 1-5 are carried out for a sufficient number of azimuth angles 0*, the resulting state 
matrices 21(00 are typically averaged to obtain the final constant state matrix A. This is 
the traditional linearization procedure, and it is also what is done in the present study for 
the calculation of the portion of the A matrix corresponding to the “average” states, i.e., for 
the submatrix A ave in equation 4.22. 

Some additional manipulations are required for the rows corresponding to the 4/rev 
derivatives x 4 p, i.e., for the submatrix A 2 1 . As for the B matrix case, first perform a 
multiblade coordinate transformation on Ar( 00, resulting in Ap(00 , and then extract the 
4/rev cosine and sine harmonics through a Fourier analysis. Define: 


± R+ = f(x R+ ,U,00 


(4.35) 



(4.36) 


A r(00 = [{Ar( 00} 1 {Ap(00} 2 • • • {Ap(00}J 


(4.37) 


A Fic = — J2 A F (00 COS 40; 

JV,/, 


i = 1 

2 N-if} 


(4.38) 


Afi s = sin 4 V7 

JV,/, 


i> i= 1 


(4.39) 
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Then it is essentially 


^21 — 


a f 4c 

a f 4s 


(4.40) 


except that the rows Ap 4c and Ap 4a must be appropriately permutated because the state 
derivative sub vector x 4 p [Eq. 4.21] is arranged with the 4/rev cosine and sine components 
interlaced rather than grouped together. 


4.2.3.2 Rows corresponding to the 4/rev states x 4 p 

The rows of the A matrix corresponding to the 4/rev state vector x 4 p, i.e., the submatrices 
A 12 and Ahhc in equation 4.22 can be obtained with the same procedure as previously 
shown for the A ave and A 21 matrices with two special treatments: 4/rev perturbation and 
kinematic relationship. 

Unlike the conventional linearization method which has a constant perturbation, the 
submatrices A 12 and A HHC are obtained by perturbing x 4 p in 4/rev frequency in both sine 
and cosine direction. The 4/rev frequency is chosen to excite the 4/rev response. Because 
the equations of motion of the helicopter are not expressed in terms of 4/rev states, they 
cannot be perturbed directly. The alternative solution is to perturb each main rotor states 
x M r of the partition x a „ e by ±Ax M r cos 4?/)* and ±Ax M r sin 4^. This is the same as 
perturbing x M r 4c and x M r 4s by ±Ax M r, respectively. 

Another important aspect regarding the 4/rev perturbation is the kinematic relationship 
between the rotor states. There are several types of kinematic relationships, and one of 
them is the integral relationship such as 


^(Ac) — Ac (4.41) 

-^(Ac) = Ac (4.42) 

Because the periodic nature of the 4/rev states, the kinematic relationships are maintained 
in a different way. Continuing with the previous example, the first and second derivatives 
of equation 4.18 with respect to time are given by: 


Ac (VO 
Ac (VO 


Ac (VO 


Ac a „ e (VO + Ac 4c (VO cos V> + Ac*. (VO sin V> (4. 1 8) repeated 
Ac a „ e + (Ac*. + 4UAc*J COS 4^ + (Ac* s - 4fiAc*.) sin 4V» 

V y > V ^ > 

^ic* c 


Pica W + (Ac*. + 8UA C4s - 16U 2 Ac*.) COS 4^ 


lc *c 


(4.43) 


+ (Ac* s - 80 A C4c - 16U 2 Ac*J sin 4-0 


l c 4 s 


(4.44) 
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Although /3± c , Pi c , and p lc on the left-hand side of equations 4.18, 4.43, and 4.44 
correspond to the integral relationships, the primed and dotted variables on right-hand side 
of equations 4.18, 4.43, and 4.44 do not maintain the same integral relationships. For 
example, 


S<A~> 

+ PL 

(4.45) 


^ PL 

(4.46) 


P 1 P'U 

(4.47) 

>U.) 

PlCAs 

(4.48) 


(4.49) 


Define all the dotted and non-dotted variables such as Pi C4c , /?i C4s , Pi C4c , and Pi C4c on the 
right-hand side of equations 4.18 and 4.43 to be the 4/rev rotor states, and the primed 
variables such as P\ CAr and P[ C4s to be the global 4/rev rotor states (for this simple example). 
Then, the kinematic relationships between the 4/rev rotor states and global 4/rev rotor states 
are as follows: 


PlCic 

— Pic4 C + 4D/3 1c4s 


(4.50) 

PLs 

= Pic ia ~ 4fZ/3i C4c 


(4.51) 

P'U 

= Lac + 8 ^/ 5 i C4s 

- 16f l?Pi CAc 

(4.52) 

PLs 

= PlC4,s — 8^l$i C4c 

- 16f2 2 /?i C4s 

(4.53) 


For this research, the 4/rev rotor state vector and the global 4/rev rotor state vector are 
shown in equations 4.21 and 4.26, respectively. 

The kinematic relationship must always be maintained throughout the linearization 
process. For instance, if the Pi c was perturbed by a constant APi c , the Pic must also 
be perturbed by ^(A/3i c ) at the same time to maintain kinematic consistency. Because the 
time derivative of a constant perturbation is zero, the traditional linearization method only 
perturbs one state at a time while the rest of the states remain fixed. 

For a 4/rev perturbation, if P lc is perturbed by APi c cos 4-0 , p tc must also be perturbed 
by A(A/3 lc cos4^) at the same time. Conversely, if Pi c is perturbed by Ap lc sin Pic 
must also be perturbed by ^(A/3i c sin 4p). It is important to remember that the equations 
of motion of the helicopter are not expressed in terms of 4/rev rotor state, x 4 p cannot be 
perturbed directly. The procedure described below perturbs each main rotor state with x M p 
in both sine and cosine direction at a 4/rev frequency. 

The calculation of submatrices A 12 and Ahhc proceeds as follow: 

For every azimuth angle -0F 

1 . Perturbation of the 4/rev cosine component 
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(a) Perturb the j th main rotor state x M n :j (V0 of partition in the state 

vector ^ave(V'i) by A s Mfij cos 40* , i.e., let the perturbed state vector x 4c+ (0*) 
be 


XMR-t 


X 4c+ (00 


XMRj - 4 + [-^tlAxMRj sin 40*] 


XMRj + Ax MR, COS 4 0* 


(4.54) 


XMR L 

where the subscript “4c+” denotes the positive 4/rev cosine perturbation in the 
central difference calculation. 


(b) Substitute the perturbed state vector x 4c+ (0*) into the system of equations 
of motion of the helicopter, to obtain the perturbed state vector derivative 
x 4c+ (0*). Because the rotor equations are formulated and implemented in the 
rotating system, x 4c+ (0*) must first be converted to the rotating system using 
a multi-blade coordinate transformation that yields the corresponding rotating 
state vector xb 4c +(0*) 


XR4c+(i) = f(xR4 c+ (0i),u(0*),0i) (4.55) 


The control vector u(0*) corresponds to the desired trim condition, and is held 
constant during the perturbation. 

If and only if x MRj is one of the displacement states (/3 0 , d ]c , 3 ]s , /3 2 , Co, Cic, 
Cis, C 2 , 00, 01c, 01 s, 02), its derivative state x MRj _ 4 (0o, $1 c, $ 2 , Co, Cic, Cis, 
C 2 , 0o, 0i c, 0is, 02) also needs to be perturbed by —AVtAxMRj sin 40* at the 
same time. This additional perturbation is represented by [. . .] in equation 4.54. 

(c) Repeat the two previous steps with a negative perturbation of the j th main rotor 
state, xmrj — Ax mu, cos 4 -0,, and build the derivative using central difference 
approximations. This derivative is the j th column of an interim matrix Pr(0i) 
at the azimuth angle 0* , that is: 

{-Pr(0*)},- ~ (x/? 4c +(0i) - xb 4c _ (0;)) (4.56) 

3 2A x MR . 


(d) Repeat the three previous steps for each of the L elements of the main rotor 
state vector partition x MR , i.e., for the main rotor state in x M r :i , j = 1, . . . , L, 
to obtain the first half of the interim matrix Pr(0O 

PrM = UPrM}, IPrM}, ■ ■ ■ IPrM}lL*l < 4 - 57 ) 
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2. Perturbation of the 4/rev sine component 


(a) Perturb the j th main rotor state x M r :i of partition x A / R (A) in the state 
vector x a „ e (^j) by A xmRj sin 4^, i.e., let the perturbed state vector x 4s+ ('0i) 
be 

x B 

XMR-L 


x 4s +(^j) 


x M Rj- 4 + [4fiA XMRj cos4-0i] 


x M Rj + a x M Rj sin 4^ 


(4.58) 


XmR l 

where the subscript “4s+” denotes the positive 4/rev sine perturbation in the 
central difference calculation. 


(b) Substitute the perturbed state vector x 4s+ (^j) into the system of equations of 
motion of the helicopter to obtain the perturbed state vector derivative x 4s+ ( 7 /^). 
Because the rotor equations are formulated and implemented in the rotating 
system, x 4s+ ('0 i ) must first be converted to the rotating system, using a multi- 
blade coordinate transformation that yields the corresponding rotating state 
vector x R 4 S+ (A) 


Xi? 4s+ (^i) = f(x fl4s+ (^), u(V>i), i/Ji) (4.59) 


The control vector u(4/y) corresponds to the desired trim condition, and is held 
constant during the perturbation. 

If and only if xmr^ is one of the displacement states, its derivative xmr :i _ a also 
needs to be perturbed by +AQAxmr.j cos 4^ at the same time. This additional 
perturbation is represented by [. . .] in equation 4.58. 

(c) Repeat the two previous steps with a negative perturbation of the j th main rotor 
state, xmrj — A xmRj sin 4^, and build the derivative using central difference 
approximations. This derivative is the (L + j) th column of the interim matrix 
P R {$i) at the azimuth angle ip l , that is: 

{PR(A)} L +j ~ TT— (x/? 4 s+(^i) - x R4s _(V\)) (4.60) 

ZlXXMRi 


(d) Repeat the three previous steps for each of the L elements of the main rotor 
state vector partition x M k, i.e., for the main rotor state in xmr :i , j = 1 , ,L, 
to complete the second half of the interim matrix Pr^) 


Pr(A) 


{PrW}i • • • {Pr(^)}l {Pr{Ml+ 1 • • • { P ^)} 2 lL 9L 

(4.61) 
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The next step is to perform a multi-blade coordinate transfoimation to obtain an interim 
matrix Pp(0j) entirely in the fixed system. Then, after steps 1-2 are carried out for one 
rotor revolution, the resulting interim matrices Pp(0) are averaged to obtain the state 
matrix A 12 . Next, the columns of A 12 must be appropriately permutated because the state 
sub vector x 4 p [Eq. 4.21] is arranged with the 4/rev cosine and sine components interlaced 
rather than grouped together. 

The state matrix Ahhc can be obtained by extracting the 4/rev cosine and sine harmonics 
from Pp(0j) using Fourier analysis. Define: 


2 ^0 


T r J2 P F (^i) cos ^i 

i = 1 

(4.62) 



T r£ i M0i) sin4 0* 

i= i 

(4.63) 


Then it is essentially 


Ahhc — 


Af 4 c 

A-f 4s 


(4.64) 


except that the rows and columns of Ap 4c and Ap 4s must be appropriately permutated 
because the state subvector x 4 p [Eq. 4.21] is arranged with the 4/rev cosine and sine 
components interlaced rather than grouped together. 

There is one last special treatment related to state matrix Ahhc • Both submatrices Ap 4c 
and Ap 4s do not contain the 4/rev state derivatives x 4 p. Recall that the interim matrix 
Pp(0i) contains the perturbed state vector derivative which has elements such as 3\ c . The 
Fourier analysis only extracts the global 4/rev rotor states (fi'( CAr and not the 4/rev 
rotor states 0i C4c or /3 4c4s ). To conform with standard state-space representation x = Ax + 
B u, the global 4/rev rotor states in Ap 4c and Ap 4g are converted to the 4/rev rotor states 
using the kinematic relationship as shown in equations 4.50-4.53. 


4.2.4 Extraction of the feedforward matrix D 

Consider that the hub loads at hub in body axes can be written in symbolic form as: 

F = g(x, x, u; 0) (4.65) 

The extraction procedure for the feedforward matrix D is the same as the one for the control 
matrix B except the subject of the interest is the hub loads F instead of the state vector 
derivatives x. 

The calculation proceeds as follows. For every azimuth angle 0,: 

1. Perturb the k - th element iikiAh) of the control vector u (pilot and HHC controls are 
treated in exactly the same way) by A u fc , i.e., let the perturbed control vector u + (00 
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be 


u+(^i) = < 


u 1 

U 2 

'U’k T All k 


(4.66) 


tt m 

where the subscript “+” denotes the positive perturbation in the central difference 
calculation. 

2. Substitute the perturbed control vector u + (0j) into the equation 4.65, to obtain the 
perturbed hub loads F + 

F + = g(x,u + ,ipi) (4.67) 

The state vector x corresponds to the desired trim condition, and is held constant 
during the perturbation. 

3. Repeat the two previous steps with a negative perturbation of the A:-th control, — 
A Uk, to obtain the perturbed acceleration vector F_. 

4. Build the derivative using central difference approximations. This derivative is the 
k - th column of the D matrix at the azimuth angle 


{£>(*)}, 


2Aul 


(F i_ — F_ 


(4.68) 


5. Repeat the four previous steps for each of the m elements of the control vector, i.e., 
for Uk, k = 1, . . . , m, to obtain the complete control matrix D(ipi) 


DM) = {{DM)}, {DM)}, ■ • • {£>(*)}„] 

6. Repeat steps 1-5 for azimuth angles , il) l for one rotor revolution. 


(4.69) 


7. Extract the average, 4/rev cosine, and 4/rev sine harmonics of D(ipi) using Fourier 
analysis. 

Define: 


D a ve 


Dac — 


Da s — 


1 


A" wo 
1=1 

(4.70) 



TT cos 4 0* 

7N V i= i 

(4.71) 

2 -^0 


TrE^W 8 ' 11 4 0* 

i= 1 

(4.72) 
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Then it is essentially 


[-D31 -D32] — D, 

[-D41 -D42] 


*4.73) 

(4.74) 


except that the rows of D± c and I)\ s must be appropriately permutated because 
the output subvector F 4 p [Eq. 4.25] is arranged with the 4/rev cosine and sine 
components interlaced rather than grouped together. 


4.2.5 Extraction of the output matrix C 

The general procedure to extract the state matrix C is similar to that of the control matrix 
A. 


4.2.5.1 Submatrix C 2 2 

The submatrix C 22 relates the 4/rev rotor states x 4 p to the global 4/rev rotor state vector 
y 4 P; i.e., C '22 is a kinematic matrix. Using 3[ cu and (3' lCis as an example, the kinematic 
equations for /3 lc are 


# C4c = /W + 40/? lc4s (4.50) repeated 

Pic 4 , s = Pic* ~ 40/3 1C4c (4.51) repeated 

(4.75) 


Re-write the above equations in matrix form: 


P'lC4c 


' 1 

0 

0 

40 ' 

. PlCis 


0 

1 

—40 

0 


/^ le 4c 
Pi C4s 
PlC4c 
Pi C4s 


Likewise, the C '22 can be structured as follows: 


Yap — 


H 0 0 

0 H 0 
0 0 H 


X4p 


(4.76) 


(4.77) 


■J0001U0 0 O' 

0/0001U0 0 

00/00 01U0 

000/0 0 0 w 


0 4U 
-4fi 0 


(4.78) 
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4.2.5.2 Submatrices C 3 1 and C 41 

The rows of the C matrix corresponding to the average states x ave , i.e., the submatrices C 31 
and C 4 1 in equation 4.23 can be obtained with the same procedure as previously shown for 
the A ave and A 2 1 matrices, i.e., through the following steps. 

For every azimuth angle -0F 


1 . Perturb the k - th element x avek (t4) of the partition x„,, e of the state vector x by Ax k , 
i.e., let the perturbed state vector x + (^j) be 


%ave i 
%ave 2 


X+(^i) 


%avek 


+ Ax k > 


(4.79) 


*£<rue jv 


X4p 


where the subscript “+” denotes the positive perturbation in the central difference 
calculation. 


2. Substitute the perturbed state vector x + (^i) into equation 4.65 to obtain the perturbed 
hub loads F + . 

F+ = g(x + , u, ipi) (4.80) 

The control vector u corresponds to the desired trim condition, and is held constant 
during the perturbation. 


3. Repeat the two previous steps with a negative perturbation of the /c-th average state, 
x avek — A Xk, to obtain the perturbed hub loads F_. 

4. Build the derivative using central difference approximations. This derivative is the 
k - th column of an interim matrix P{^i) at the azimuth angle i ) % , that is: 

{m)} * “ 2^ (F+ “ F - ) <4 ' 81) 


5. Repeat the four previous steps for each of the N elements of the state vector partition 
x aiJe , i.e., for x avek , k — 1, . . . , N, to obtain the complete matrix P(^i) 


jv] (4.82) 


6. Repeat steps 1-5 for azimuth angles for one rotor revolution. 
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7. Extract the average, 4/rev cosine, and 4/rev sine harmonics of using Fourier 
analysis. 

Define: 


Then it is essentially 


r 

w ave 

C 4c 

c 4s 


1 ^0 


F EPW 

n * i,l 

(4.83) 



TrE P W c °s 4 i 

i= i 

(4.84) 



T ^^ P (^) sin4 ^i 
i=1 

(4.85) 


C 31 = c, 

cv, = 


'ave 

C 4c 

c 4s 


(4.86) 

(4.87) 


except that the rows C 4c and C 4s must be appropriately permutated because the output 
subvector F 4 p [Eq. 4.25] is arranged with the 4/rev cosine and sine components 
interlaced rather than grouped together. 


4.2.5.3 Submatrices C 32 and C 42 

The submatrices C : > 2 and C 42 can be obtained with the same procedure as previously shown 
for the A 12 and A HHC matrices which proceeds as follow: 

For every azimuth angle -0d 

1 . Perturbation of the 4/rev cosine component 

(a) Perturb the j th main rotor state x M n :/ (A J i) of partition x MR (?/y) in the state 
vector x a „ e (-0i) by Ax M r :i cos 4 fa, i.e., let the perturbed state vector x 4c+ (^) 
be 


%MR i 


X 4c +('0i) 


xmr, 4 + [-4DA x M r, sill 4^] 

< . 


xmr, + a x M r 3 COS 4 i>i 


(4.88) 


XmRl 

where the subscript “4c+” denotes the positive 4/rev cosine perturbation in the 
central difference calculation. 
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(b) Substitute the perturbed state vector x 4c+ (0j) into equation 4.65 to obtain the 
perturbed hub loads F 4c+ (0j). 

F 4c +('0i) = g(x 4c+ (0j), u (0*), ipi) (4.89) 

The control vector u(0j) corresponds to the desired trim condition, and is held 
constant during the perturbation. 

If and only if xmr^ is one of the displacement states (J3 0 , 3\ c , 6\ s , f3 2 , Co, Cic, 
Cl-. C2, 00. 01 c. 0i-, 02), its derivative state x M r 3 _ a (P 0 , Pi « Pis, P 2 , Co. Ci« Ci-, 
C2, 0o, 0i c, 0is, 02) also needs to be perturbed by — 40A sin 4i/y at the 
same time. This additional perturbation is represented by [. . .] in equation 4.88. 

(c) Repeat the two previous steps with a negative perturbation of the j th main rotor 
state, xmRj — AxmRj cos4i/y, and build the derivative using central difference 
approximations. This derivative is the j th column of an interim matrix P('iPi) at 
the azimuth angle , that is: 

(F 4c+ (0O - F 4c _(0,)) (4-90) 

(d) Repeat the three previous steps for each of the L elements of the main rotor 
state vector partition x M k, i.e., for the main rotor state in xmr : i , j = 1 ,L, 
to obtain the first half of the interim matrix P(ipi) 

P(’Pi) = l{P(’Pi)},{P(’Pl)}2---{P(’Pi)} L l xL < 4 '«) 


2. Perturbation of the 4/rev sine component 


(a) Perturb the j th main rotor state x M r :i i'Ct) of partition x MK (0j) in the state 
vector n.ave (0i ) by A xmrj sin 4 0*, i.e., let the perturbed state vector x 4s+ (0j) 
be 

%MRx 


X 4s+ (0j) = < 


x M R, 4 + [4fiAaiMi?,. cos 40j] 


+ Aimr, sin 40* 


(4.92) 


%MRl 

where the subscript “4s+” denotes the positive 4/rev sine perturbation in the 
central difference calculation. 
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(b) Substitute the perturbed state vector x 4 . s+ (?/;,) into equation 4.65 to obtain the 
perturbed hub loads F 4s+ (^j). 

F 4s +('0i) = g(x 4 s+(^), u(V>i), A) ( 4 -93) 

The control vector u(-0j) corresponds to the desired trim condition, and is held 
constant during the perturbation. 

If and only if xmrj is one of the displacement states, its derivative xmr j _ 4 also 
needs to be perturbed by +40 A xmrj cos 4 fa at the same time. This additional 
perturbation is represented by [. . .] in equation 4.92. 

(c) Repeat the two previous steps with a negative perturbation of the j th main rotor 

state, xmr 3 — sin 4)+, and build the derivative using central difference 

approximations. This derivative is the ( L + j) th column of the interim matrix 
P(ipi) at the azimuth angle ?+ , that is: 

{m)} i+ ; « 2^- (F ls+ (*) - F 4 ,_(*)) (4.94) 


(d) Repeat the three previous steps for each of the L elements of the main rotor 
state vector partition x M k, i.e., for the main rotor state in xmr : i , j = 1 ,L, 
to complete the second half of the interim matrix P(ij)i) 


P(A) 


{P(1>i)}i ■ ■ • {P^i)} L {P(^)} L+ i ■ ■ ■ {p(A)} 2L 


(4.95) 


3. Repeat steps 1-2 for azimuth angles ?+ for one rotor revolution. 

4. Extract the average, 4/rev cosine, and 4/rev sine harmonics of P{^i) using Fourier 
analysis. 

Define: 


Then it is essentially 


C a ve 

= 

A ib N * 
iV V' i=i 


(4.96) 

c, c 

= 

c\ 

TF E p tyi) 

i=i 

0 

0 

Cfi 

(4.97) 

Cas 

= 

c\ Nijj 

N * i.i 

sin 4^j 

(4.98) 


C32 

r 

w aue 


(4.99) 


C42 

c 4c 

C 4s 


(4.100) 
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except that the rows and columns of C 4c and C' 4 . s must be appropriately permutated 
because the state subvector x 4 p [Eq. 4.21] and the output subvector F 4 p [Eq. 4.25] 
are arranged with the 4/rev cosine and sine components interlaced rather than 
grouped together. 

4.3 Application to simple rotor equations 

In this section, the perturbation technique described in the previous section is applied to 
a simple example, namely the flap equation of motion of a 4-bladed isolated rotor written 
in fixed-system coordinates. The blades are assumed to be rigid and hinged at the axis of 
rotation. This simplified model represents a useful test case because states and harmonics 
appear explicitly in the equations of motion, and therefore can be manipulated directly, 
rather than being hidden in the more complicated numerics of the model used in the 
remainder of this research. 

The flapping equations of motion in the rotating system for a 4-bladed rotor with rigid 
blades hinged on the axis of rotation, and flapping degrees of freedom only can be expressed 
as 


A + A A = 


7 -^cosA Q + ^sinA^ A - ( 

Q + ^ sin A + ^/i 2 sin 2 a) 0* 

/i sin A^) A 


8 6 


/i sin A A 


(4.101) 


where A, 66 are the main rotor inflow and the blade pitch angle, respectively. After 
performing the multi-blade coordinate transformation, the equations of motion are: 


where 


A) 


A 


A 


’ A ‘ 

Ac 

+ 67 

Ac 

+ A" 

Ac 


A 

A s 

As 

As 


A 

. A . 


. A . 


. A . 


. ^4 . 


C 


l 

8 

0 

6^ 

0 



A/i sin 2?/; 


Tit 1 

2 


A+ cos 2A 


0 

— sin 2 A 

^/i cos 2A 

2 . 

8 


(4.102) 


(4.103) 
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I< = 


IV 

0 


. —'Ip 2 sin 20 


— 1 + v 2 + ^p 2 sin 40 
-| - ^P 2 cos 4 ^ + iqP 2 
— jP cos 2 0 


0 

| - i \p 2 cos 4-0 + jqP 2 
— 1 + u 2 — jLp 2 sin 4 0 
— | p sin 20 


-|/i 2 sin 2-0 
-|/icos20 
- 2 // sin 20 
v 2 

(4.104) 


F\ = 


7 


Jr /As 

6 

7 




0 + 


77 (l + /l 2 ) 6*4 C + -t-pOss — ^pO 3, 


7 

6' 


cos 4-0 


^ (l + / i2 ) $4s + -^p^3c ~ 


sin 40 


(4.105) 


Fo = 


+ 


? (l + 7}P 2 ) °lc - T^P 2 ® 


16 ' 


3c 


7 .. 2 /, , 7 , 1 . , 7 ^ 1.2 


6\ c + —(i + -p )e, c + -[i + -^ 


~>3c 


1 ( l + ^ 2 ) ^ | 


1 + 2 ^ ) ® 3s 


7 2 . ' 
16^ 9ls . 


cos 40 
sin 40 


1 p 2 9 5s sin 80 — -/-p 2 9 5c cos 80 


16 


16 


(4.106) 


F» = 


4 A/i + 8 ( X + ° ls + 3^° ~~ \6^ 0zs 


7 


3 

2' 


77 ( 1 + 77 P~ ) $3s + % ( 1 + ~ p 2 ^\ 9§ s + -J-p 2 9\ s + 0/C0 4c 


3 

2 1 


16 ' 


7 

3 ; 


^/i6* 4s — ^/ f2 ^ic — g (■*■ + 2 ^ 2 )^ 5c ■*" g (1 + 2^ 2 ) ^ 3c 

p 2 9 5c sin 80 H p 2 9 5s cos 80 


16 


16 


cos 40 
sin 40 

(4.107) 


F 4 — (——p 93s + 7:p9is + t:P~ 9 q + tt;P 2 9 ac \ cos 20 

V 6 6 8 16 / 

+ _ ^A^ic + e /u6 ' 3c ) sin 

+ f^^ 4 * - sin 60 + (y^/^c + ^r/^Ss) cos 60 ( 4 . 108 ) 
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To illustrate the technique for the extraction of a linear model that included 4/rev 
characteristics, the longitudinal flapping equation of motion in the fixed system is: 


Pic + 2/3i s + (u 2 — l) Pic 


-j-l? sin Ail) Pic + J/i sin 2 00 2 - J Pic 

16 6 6 8 

Q/i cos 2i/?j p 2 - Q - cos 4 0 + Pis 

(i + ^ 2 ) 0 


7 


1 c r 1^ ^3c 

lo 


7 , ,2/ 


7 


16 


— — A* 01c + — 1 + — At 05c + — 1 + —ft 0, 


7 


7 3c 


7 


7 


1 + 70 05s + 77 ( 1 + 70 03s — T0 0 


7 


16 


7 ls 


cos 4-0 
sin 4-0 


- 70 2 0 5s sin 80 - -^ 0 2 0 5c cos 80 
16 16 


(4.109) 


For simplicity, only the — ^/r 2 sin 40/0 term from aerodynamics and |(l+^)(0 3c cos 40+ 
03 S sin 40) terms from the 3/rev input, and — ^0 2 0ts sin 40 term from the longitudinal 
cyclic pitch angle are retained. The rest of the variables are represented by M e and Mg 
terms. Equation 4.109 can therefore be rewritten as: 

2 

Pic = - ~0ic + ( V 2 - sin 40 - l) /0 

+ y ( X + ^ /i2 ) (^ 3c cos + 6,35 sin 4 ^) _ 16 /i2 ^ ls sin + M e + Mg 

(4.110) 


4.3.1 Prescribed solution form 


The assumed solution for equation 4.110 has an average plus 4/rev cosine and sine 
components as shown in equation 4.18. Substituting equations 4.18 and 4.43 into 
equation 4.1 10, the longitudinal flapping equation of motion becomes: 


Pic — — 


7 


Picave + (0ic 4c + 4f/0ic 4s ) cos 40 + (/0 4s - 4O/0 4c ) sin 40 


+ ( 


1/2 _ sin4 ^ ~~ ^ (/ lcw + ^ lc4c cos4 ^ + ^ lc4s sin4 

,2 


+ - (l + y) i?3c cos 40 + e 3s sin 40) - :^0 2 0 i s sin 40 + M e + Mg 


(4.111) 


Note that equations 4.43 and 4.44 contain an average value and harmonics at only 4/rev, 
which results from the original assumed solution defined in equation 4.18. However, the 
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equation of motion shown in equation 4.111 also has 8/rev frequency components that 
result from the aerodynamic term in equation 4. 1 10. 


J£_ 

’ 16 


sin Aif) (3 lc = 


7 A* 

’ 16 
11? 
’ 16 


sin Aip(f3 1Cave + f3 1C4c cos 4 ^ + f3 lC4s sin 4 


(Picav* sin 4^ + ^ + \pic ic sin 8 ^ - \pic±s cos 8 


(4.112) 


Therefore, equation 4.111 can also be written as follows: 


Pic — 


PlCave + (/3lc 4c + 40^1^) COS 4^ + (/3 ic4s - 4 ^/3ic 4c ) sin 4^» 

(z / 2 - l) (/5i Ca „ e + /3 i C4c cos 4-0 + /3 1c4s sin Aipj 


7p 
16 
7/i 2 /1 


(/?i<w sin 4'0 + ^/3 1c4 .,) 


16 (g/W sin 8 V> - COS 8 ^) 

2 


^ (1 + —) (& 3 C COS 4^ + d 3s sin 4 

-p '~H 2 0is sin 4^ + M e + M 0 
16 


(4.113) 


4.3.2 Perturbation of the equations of motion 

To simplify the expression, the perturbations of the Mg and Mp terms in equation 4.1 1 1 are 
not shown here, but are not eliminated from equation 4. 1 1 1. 


4.3.2.1 Perturbing /3 lc by ± A/3 lc at I/?* 


Putyi) 


01c 


7 


PlCave + Pic 4c + 4 ^lc 4s COS Alp 


+ 


2 

Pic*. ~ 4 ^A c ic sin 4^) | + (v 2 - sin Aip - 1 


x 

+ 


(/?icw ± A/3i c ) + /3 i C4c cos 4^ + Pic4 S sin 4^ 
2 

|(y + l) ( 6*30 cos 4^ + 6 > 3s sin 4^) 

^~li 2 Q\ s sin 4^ 

16 


(4.114) 


where the superscript ± represents the direction of the perturbation, and the subscript Pi Cave 
represents the perturbed state variable. 
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4.3.2.2 Perturbing (3 lc by ± A (3 lc cos 4-0 at 0* 


01c(0i 


± 

/3ic 4c 


^\Pl Cave + (/W + 4 ^/3lc 4s ) COS 40 


+ 


X 


2 

0i C4s - 4fi(/3i C4c ± A0 ic ) sin 40 J + (v 2 - ~^r sin 40 
Pi c a „ e + (/3 i C4c ± A0 lc ) cos 40 + 0i C4s sin 40 


16 


+ |(y + l)(6>3 C cos40 + 0 3s sin4 


- -^/x 2 6>i s sin 4 0 

16 


4.3.2.3 Perturbing f3 lc by ± A f3 lc sin 40 at -0* 


A c (^i 


± 

^ lc 4s 


+ 


X 


7T *{ /^ltW + 


Pi C4c + 4 ^ (Ac 4s ± A/3l c ) 


COS 4 0 


Pu 4s - 40/3, 


^4c 


sin 40 j + ( 


2 7/^ • 

z/ sin 

16 


4-0 - l) 


Cleave + /3 i C4c cos 40 + (0 1C4s ± A0 lc ) sin 40 


+ ^ ( ■ y + 1 ) (6*3c cos 40 + 6» 3s sin 40) 


O0-2 
16 


fi 0i s sin 40 


4.3.2.4 Perturbing /3 lc by ± A/3i c at 0* 


$lc(0i) 


± 

/^Icane 


7 


{PlCave ± A01 C ) + 


Pic 4c + 41101, 


^4s 


cos 40 


+ 


X 


0ic 4s — 4H0 lc 


sin 


2 

40 j + (0 - sin 40 - l) 


PlCave + PlC4,c COS 4 V> + Plc is SU1 40 
2 

+ |(y + l) (^ 3 c cos 40 + 0 3s sin 4 

- -^/i 2 6»i s sin40 

16 


(4.115) 


(4.116) 


(4.117) 
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4.3.2.5 Perturbing /3 lc by ± A/3i c cos 4-0 at 0* 


hctyi) 


± 

@ lc 4c 


+ 


X 


7T /^lCave + 


/3ic 4s -4^A 


(pic4c =t A/?i c ^ + 40/3i 
in 40 j + ( 


2 7/^ • 

z/ sin 


cos 40 
4-0 — l) 


Cleave + /5lC4c COS 4 0 + Plc is 8W. 40 
2 

+ |(^ + l)(0 3 cCos40 + 0 3s sin40) 


_ 7 _ ! ,2 
16 


H 9 \s sin 4-0 


4.3.2.6 Perturbing /3 lc by ± A/3i c sin 40 at 0* 



r> | /^1 Cave ~ 

0i C4c + 4f20i C4s 


frc 4s 8 l 

L -1 


+ 


x 


IC 4 c 


(0ic 4s ± A 0 lc ) - 4 Q /0 
Ae.™ + Ac 4c COS 40 + Ac*. sin 4 0 


cos 4 0 

2 

in 4-0 1 + (V 2 — y- sin 40 


J -Cave 

,2 


+ |(y + l) ( 0 3c cos 40 + 0 3s sin 40 ) 


7^At 2 0i s sin 40 
16 


4.3.2.7 Perturbing 0 is by ± A0 is at 0* 


£l c (0i) 


± 

O I fi^-Cave 

9is 

8 l L 


01 C 4c + 41201, 


C4s 


cos 40 


+ 


X 


filers — 4110! 


C4c 


sin 40 1 + ( 


2 7 V ■ A I ! 

Z2 tt- sin 40 — 1 


16 


k Cave + Ac 4c cos 40 + 0 iC4s SHI 40 
2 

0 3c cos 40 + 9 3s sin 40 


^ + l) 


8 V 2 

- y A < 2 (^ls ± A 0 ls ) sin 40 


(4.118) 


(4.119) 


(4.120) 
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4.3.2.8 Perturbing 0 3c by ± A 0 3c at 0* 


Plctyi) 


+ 


X 


7H Aca« e + 


/3ic 4c + 4f2/3i 


^4s 


cos 4 0 


(3 1C4s -M(h 


^4c 


sm 


2 

40 j + (z/ 2 — y^- sin 40 — 1 


Pic ave + Ac4 C cos 4 0 + Ac*. sin 4 0 
2 r 

|(y + l) ( 6*30 ± A 6 » 3c ) cos 4 0 + 6 > 3s sin 40 

-^/i 2 6 »i s sin 4-0 
16 


(4.121) 


4.3.2.9 Perturbing 0 3s by ± A0 3s at 0* 


Plcfyi) 


03s 


7 


Pi Cave + 


/3ic 4c + 41201, 


^4s 


cos 40 


+ 


X 


0ic 4s - 4001 


C 4c 


sin 40 1 + ^ 


2 7 /* • A , . 

1/ t— Sill 40—1 


16 


Acave + PlC4o COS 4 ^ + Plc 4 s Sil1 4 0 


|(y + l) # 3c cos 40 + (0 3s ± A 6 » 3s ) sin 40 

TZl^Qis sin 40 
16 


(4.122) 


4.3.3 Extract four/rev harmonic components 

The average value and the 4/rev harmonic components of the perturbed equations 4.114- 
4.122 can be obtained by applying the Fourier series approximation over one sample 
window. The length of the sample window used in this study is one rotor revolution or 

0 = 0-27T. 
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The elements of the state matrix can be calculated using central difference approximation. 
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4.3.4 Construction of the state equation 


Combining equations 4.123-4.149 and equations 4.50-4.51, the state equation of the 4/rev 
LTI-HHC model can be constructed as follows: 
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The vector on the left-hand side of equation 4.150 consists of both primed and dotted 
variables. To conform with the standard state-space representation, x = Ax + Bu, 
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the primed variables are replaced with the dotted variables using equations 4.50-4.53 as 
follows: 
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Substitute equation 4.151 in equation 4.150, and re-arrange the equation. The 4/rev LTI- 
HHC model is given by: 
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4.3.5 Analytical model validation 

From equation 4.152, Pi Cave , Pi C4c , and (3 lCig can be expressed as follows: 
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Substitute equations 4.155-4.157 in equation 4.44, 
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Compare equation 4.158 with the equation 4.1 13 which is shown below: 
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(4.113) repeated 


The differences are the 8/rev frequency contents which were truncated from the 
linearization procedures. These 8/rev frequency contents can be retained in the linear model 
if the prescribed solution of equation 4.18 also contains 8/rev components as well as 4/rev 
components. 
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4.4 LTI-HHC model validation 


In this section, the LTI-HHC model was validated against the full-blown nonlinear 
helicopter model by comparing their 4/rev hub load and the rotor responses over several 
flight configurations. 

4.4.1 4/rev hub load comparison 

Validation was conducted for forward velocity of 40, 80, and 120 kts. Each case starts 
from the trim condition without HHC input. After two rotor revolutions, the HHC input is 
engaged, and results are shown in figures 4.1-4.15. The HHC input is a 3/rev input with an 
amplitude of 0.6° at 0° phase angle. 

The 4/rev hub loads calculated from the nonlinear helicopter model are the output of 
the harmonic analyzer which contain time delays. On the other hand, the 4/rev hub loads 
calculated from LTI-HHC model are obtained instantaneously. There is no time delay 
associated with the sample window as with the harmonic analyzer. The effects of the 
sample window can be approximated by an equivalent lowpass filter that must be included 
in the output of the LTI-HHC model before it is compared to the nonlinear 4/rev results. 

As illustrated in the figures 4.1-4.15, the LTI-HHC model produces the levels of 4/rev 
vibrations that are very close to the nonlinear 4/rev vibration. A close match is seen not 
only in steady state condition but in transients. The small ripples in the nonlinear results 
are the 8/rev and higher frequencies that were not modeled in the LTI-HHC model. 

The effect of pilot input on vibrations is illustrated in figures 4.16-4.30. The input is 
a lateral cyclic doublet input with an amplitude of one stick inch. Lrom the figures, the 
LTI-HHC model shows the capability of predicting the 4/rev hub loads. There are strong 
8/rev and higher frequencies in the nonlinear results that are not modeled in the LTI-HHC 
model. However, the 8/rev frequency can be captured by the LTI-HHC model if the 8/rev 
frequency is prescribed in the assumed solution. The dimension of the LTI-HHC model 
matrices will increase to accommodate the additional 8/rev rotor states. 

4.4.2 Rotor states comparison 

Ligures 4.31-4.33 compare the rotor states from both the nonlinear and the LTI-HHC model 
simulation for the 120 kts case. The HHC input is the same 3/rev input in the previous case. 
The rotor states compared in the figures are full values, i.e., they are not 4/rev rotor states. 
The rotor state data of the nonlinear helicopter model is obtained directly from the time 
integration of the equations of motion. Since this set of data has never passed through the 
harmonic analyzer to obtain its 4/rev components, the effect of the sample window is not 
included. 

To compare with nonlinear results, rotor state data of the LTI-HHC model is constructed 
by modulating instantaneous 4/rev rotor state data as shown in equation 4.18 without 
including the effect of the sample window. These figures illustrate that the prediction from 
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the LTI-HHC model is very similar to the nonlinear helicopter model in both steady-state 
and transient condition. 

Figures 4.34-4.36 show the effect of lateral pilot input on blade rigid flap, rigid lag, 
and torsion modes for both nonlinear and the LTI-HHC models. Input is a lateral cyclic 
doublet input with the amplitude of one stick inch. Variation of the 4/rev rigid flap and 
4/rev rigid lag modes within the nonlinear helicopter model are relatively small compared 
with variation of their mean value. With the LTI-HHC model, 4/rev characteristic of the 
torsion mode of nonlinear helicopter model as shown in figure 4.36 is predicted not only in 
the steady-state condition but also in the transient. 
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Figure 4.1. Longitudinal hub shear comparison; V=40 kts, W=14,000 lb, A 3 = 0.6°, 
03 = 0°. 
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Figure 4.2. Lateral hub shear comparison; V=40 kts, W=14,000 lb, A 3 = 0.6°, 0 3 = 0°. 
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Figure 4.3. Vertical hub shear comparison; V=40 kts, W=14,000 lb, A 3 = 0.6°, </> 3 = 0°. 
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Figure 4.4. Longitudinal hub moment comparison; V=40 kts, W=14,000 lb, A 3 = 0.6°, 

03 = 0°. 
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Figure 4.5. Lateral hub moment comparison; V=40 kts, W=14,000 lb, A 3 = 0.6°, 0 3 = 0°. 
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Figure 4.6. Longitudinal hub shear comparison; V=80 kts, W=14,000 lb, A 3 = 0.6°, 
03 = 0°. 
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Figure 4.7. Lateral hub shear comparison; V=80 kts, W=14,000 lb, A 3 = 0.6°, </> 3 = 0°. 
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Figure 4.8. Vertical hub shear comparison; V=80 kts, W=14,000 lb, A 3 = 0.6°, (p 3 = 0°. 
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Figure 4.9. Longitudinal hub moment comparison; V=80 kts, W=14,000 lb, A 3 = 0.6°, 

03 = 0°. 
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Figure 4.10. Lateral hub moment comparison; V=80 kts, W=14,000 lb, A 3 = 0.6°, 

03 = 0°. 
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Figure 4.11. Longitudinal hub shear comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 
03 = 0°. 
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Figure 4.12. Lateral hub shear comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 0 3 = 0°. 
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Figure 4.13. Vertical hub shear comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 0 3 = 0°. 
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Figure 4.14. Longitudinal hub moment comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 

03 = 0°. 
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Figure 4.15. Lateral hub moment comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 

03 = 0°. 
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Figure 4.16. Longitudinal hub shear comparison; V=40 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.17. Lateral hub shear comparison; V=40 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.18. Vertical hub shear comparison; V=40 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.19. Longitudinal hub moment comparison; V=40 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.20. Lateral hub moment comparison; V=40 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.21. Longitudinal hub shear comparison; V=80 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.22. Lateral hub shear comparison; V=80 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.23. Vertical hub shear comparison; V=80 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.24. Longitudinal hub moment comparison; V=80 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.25. Lateral hub moment comparison; V=80 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.26. Longitudinal hub shear comparison; V=120 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.27. Lateral hub shear comparison; V=120 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.28. Vertical hub shear comparison; V=120 kts, W=14,000 lb, 1-inch lateral cyclic 
doublet input. 
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Figure 4.29. Longitudinal hub moment comparison; V=120 kts, W=14,000 lb, 1-inch 
lateral cyclic doublet input. 
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Figure 4.30. Lateral hub moment comparison; V=120 kts, W=14,000 lb, 1-inch lateral 
cyclic doublet input. 
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Figure 4.31. (3 comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, </> 3 = 0°. 
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Figure 4.32. ( comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, 0 3 = 0°. 
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Figure 4.33. </> comparison; V=120 kts, W=14,000 lb, A 3 = 0.6°, </> 3 = 0°. 
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Figure 4.34. (3 comparison; V=120 kts, W=14,000 lb, 1-inch lateral cyclic doublet input. 


113 



Cic 

(deg) 



Cis 

(deg) 



Figure 4.35. ( comparison; V=120 kts, W=14,000 lb, 1-inch lateral cyclic doublet input. 
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Figure 4.36. (p comparison; V=120 kts, W=14,000 lb, 1-inch lateral cyclic doublet input. 
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5 HHC and AFCS Interaction Study 


A linear time-invariant state-space approximation that accurately models the coupled rotor- 
fuselage dynamics, including the higher harmonic response of the rotor, has been developed 
in chapter 4. This work allows several important questions to be answered regarding 
the dynamic interaction between Automatic Flight Control System (AFCS) and High 
Harmonic Control (HHC), including the effect on handling-qualities. The key breakthrough 
is in the method to extract a linear time-invariant model that includes a harmonic analyzer 
and allows the periodicity of the helicopter response to be captured. The coupled high- 
order linear model provides the needed level of dynamic fidelity to permit study of AFCS 
and HHC interaction. 


5.1 Effect of a fixed HHC input on rigid body dynamics 

To understand the potential coupling between AFCS and HHC, an analysis was first 
performed in the open-loop system to determine whether a fixed HHC input had any direct 
effect on the rigid-body dynamics. Any influence from the HHC will be indicated by the 
changes in the frequency response. Before proceeding with any interaction analyses, it is 
important to validate the baseline (HHC-off) cases of both LTI-HHC and nonlinear models 
by comparing their frequency responses against the flight test data. 


5.1.1 Open-loop frequency response validation 

In section 4.4, the LTI-HHC model was validated against the nonlinear model by comparing 
the hub load responses over several flight configurations. It is a time domain comparison, 
and it is sufficient for checking the aeromechanic quantities. For flight dynamics analysis, 
it is more common to perform the comparison in frequency domain. Figure 5.1 shows the 
P/ 5i a t frequency response comparison between the LTI model, the nonlinear model, and 
the flight test. Unless noted otherwise, all the results presented in this chapter have the 
weight of 14,000 lb at a speed of 120 kts. The frequency response of the nonlinear model 
was obtained by performing frequency sweeps in pilot lateral stick input and recording the 
vehicle roll rate response time history. The P/ 8i a t frequency response was identified by 
extracting the information from the time history data using CIFER® (ref. 58). 

Since the LTI-HHC model is already in the linear system, its frequency response can 
be calculated directly from the LTI-HHC model. Figure 5.1 shows that all three cases 
agreed with each other in the frequency range of 2-20 rad/sec. There were some small 
disagreements in the frequency range of 1-2 rad/sec between the flight test result and the 
analytical results, but the difference is not significant. Comparing the nonlinear and LTI- 
HHC frequency responses, there is also a little difference, and most of the difference is in 
the phase curve below 2 rad/sec. 
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5.1.2 Effect of an optimum three/rev input on rigid body dynamics 

Figure 5.2 shows the effect of a fixed HHC input on the rigid body dynamics for the 
nonlinear model. The fixed HHC input chosen is an optimum 3/rev input which is 
calculated from the optimization procedure that minimized the norm of 4/rev in-plane hub 
shears. The optimization procedure is similar to the one described in section 2.8. This 
figure indicates that the optimum 3/rev input has no effect on the rigid body dynamics in 
the frequency range of interest. Figure 5.3 shows the same conclusion for the LTI-HHC 
model. 

The frequency response of the nonlinear model with the optimum 3/rev input is extracted 
using the same method as the nonlinear baseline (HHC-off) case stated earlier. For the LTI- 
HHC model, one cannot simply include an optimum 3/rev input and compute the frequency 
response because the linear model will only respond at the same frequency as the input 
signal. In this case, the input signal is a 3/rev (81 rad/sec for UH-60) and it is beyond the 
frequency range of interest. To see the effect of the optimum 3/rev input on rigid body 
dynamics, one must engage the HHC loops and let the effects of the 3/rev input propagate 
through the HHC feedback loops. 

Although the results above show that the HHC input has no effect on rigid body dynamics 
(or AFCS), it does not necessarily mean the AFCS has no effect on the HHC. There is still 
a possibility that the AFCS affects vehicle vibration and indirectly affects the HHC. This 
closed-loop analysis is discussed in the next section. 

5.2 Interaction of HHC and AFCS 

A SIMULINK® simulation of the combined flight and higher harmonic control system 
was developed for analysis and optimization in the Control Designers Unified Interface, 
CONDUIT® (ref. 59). The key elements of the simulation are illustrated in figure 5.4, 
and they are: 

1. Higher-order linear airframe model that provides the flight mechanics and 4/rev 
vibration responses to both pilot and HHC inputs. 

2. Automatic Flight Control System loops based on a simple proportional-integral- 
derivative (PID) controller in roll, pitch, and yaw. 

3. Typical actuator/sensor filter dynamics. 

4. Equivalent harmonic analyzer approximates the sample window dynamics and 
equivalent time delay. 

5. Higher harmonic controller based on fixed T-matrix feedback. 

6. Zero-order-hold approximation simulates the discrete HHC update time delay. 
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Like the open-loop analysis, it is important to validate the closed-loop model to ensure 
that the linear continuous time domain model implemented in SIMULINK® is equivalent 
to the nonlinear multi-rate model. This can be accomplished by comparing the broken 
control loop response of both models. 

5.2.1 Broken control loop response validation 

Figure 5.5 illustrates the schematics of both linear and nonlinear simulation models. The 
simulation model shown in figure 5.5a is a nonlinear multi-rate system. Because the control 
system analysis was performed in the linear continuous-time domain, the entire nonlinear 
multi-rate system was converted to an equivalent linear continuous time system as shown 
in figure 5.5b. The harmonic analyzer is now embedded within the LTI-HHC model. The 
effect of the sample window is modeled by an equivalent lowpass filter. The discrete HHC 
controller is transformed to a continuous-time domain HHC controller. The discrete HHC 
update (zero-order-hold) is approximated by a Pade function. 

The broken control loop response is a method of studying loop stability; it allows one to 
determine the gain and phase stability margins. Using 0 3c broken control loop for instance, 
it is the d 3c response at point B in figure 5.5 with respect to the d 3c input at point A while 
the 3/rev-cosine and 3/rev-sine loops are open. The flight control system is also disabled 
during the frequency sweep. For the purpose of the validation, six broken control loop 
responses (3/rev-cosine, 3/rev-sine, 4/rev-cosine, 4/rev-sine, 5/rev-cosine, 5/rev-sine) were 
extracted from each model, and the direct comparisons are shown in figures 5. 6-5. 8. In 
these figures, the frequency response of the LTI-HHC model matches very well with the 
one from the nonlinear model in both the magnitude and phase curves for all six loops 
within frequency range of interest. This indicates that the linear continuous time domain 
model in figure 5.5b is equivalent to the nonlinear multi-rate model in figure 5.5a. 

Although HHC input operates at 3, 4, 5/rev frequencies (or 81, 108, 135 rad/sec for 
UH-60 helicopter), the crossover frequency of each HHC loop is only about 1 rad/sec. The 
crossover frequency, gain margin, and phase margin of each HHC loop are tabulated in 
table 5.1. Because of the high HHC input frequency, one would expect a large frequency 
separation between the flight control and HHC system and assume these two systems do 
not interfere with each other. However, research results show that not only do the HHC 
loops operate at a much lower frequency, but they are also within the frequency range of 
the flight control system. This is another indication of potential HHC/AFCS interaction. 

5.2.2 Optimization of AFCS (HHC-off) 

One way to see whether the closed-loop HHC system has any effect on the AFCS or 
handling-qualities is to optimize the AFCS for the satisfactory (Level 1) handling-qualities 
with the HHC loops disengaged (fig. 5.9). Any influence introduced by closing the HHC 
loops will be indicated by the change in handling-qualities. The AFCS implemented in 
this study is based on a simple PID controller (fig. 5.10) in roll, pitch, and yaw axis. The 
PID controller computes individual actuator command with respect to the changes in rigid 
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body states and pilot inputs. The actuator is a second order model (fig. 5.11) including 
both the position and rate saturation limits. The actuator design parameters are tabulated in 
table 5.2. 

First, CONDUIT® was used to optimize the PID gains of the AFCS, with the HHC 
loops disengaged. The PID gains were tuned to achieve satisfactory handling-qualities, 
based on the Aeronautical Design Standard (ADS-33E (ref. 60)), and standard control- 
system design specifications list below (Appendix A): 

• Eigenvalue real part (EigLcGl) 

• Crossover frequency (CrslnGl) 

• Stability margins (StbMgGl) 

• Bandwidth (BnwRoF3) 

• Step response damping ratio (OvsAtHl) 

• Crossover frequency (CrsMnGl) 

• Eigenvalue damping ratio (EigDgGl) 

• Step response rise time (RisTmGl) 

CONDUIT rapidly tuned the PID gains to achieve satisfactory (Level 1) requirements 
with minimum over-design as shown in figure 5.12. The optimized PID gains are tabulated 
in table 5.3. Each symbol in figure 5.12 represents the result for a particular loop and 
shows that all the responses lie in the light region (Level 1). For example, note that the 
roll bandwidth is 3 rad/sec [fig. 5.12d] which meets ADS-33E. The PID gains of the roll 
and yaw loops yield bandwidths in excess of the requirement in order to meet some of the 
other specifications. It is important to mention that this set of PID gains is not the best from 
the handling-qualities point of view. It is simply the lowest gains needed to satisfy all the 
design specifications while staying in the level- 1 region. 

5.2.3 Nominal T-matrix controller 

Next, the T-matrix HHC loops were engaged with a nominal gain of k= 1 (same in all six 
loops) as shown in figure 5.13. This is referred to as the “nominal” case. With both AFCS 
and HHC loops closed, the CONDUIT® HQ design specifications were re-evaluated 
without changing the PID gains. The results are presented in figure 5.14 which shows 
that the closing of the HHC loops had a negligible effect on the AFCS performance and 
overall handling-qualities. This indicates the lack of dynamic coupling of HHC into flight 
control. Therefore, no re-tuning of the AFCS was needed for the combined AFCS/HHC 
system. The lack of interaction from HHC to AFCS is consistent with the earlier system 
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identification results obtained in section 5.1.2, which also showed no effect of an HHC 
input on the rigid-body dynamic response. 

In terms of suppressing the steady state vibration level, the nominal T -matrix controller 
can reduce the vibration by a large amount. Figures 5.15-5.18 show the changes in the 
lateral and longitudinal 4/rev vibration level with respect to the HHC and pilot stick inputs 
for both the baseline case (HHC-off) and nominal (k= 1) case. At t=0, the vehicle starts 
from a steady state condition, and the 4/rev vibrations are maintained at a steady level. The 
baseline 4/rev vibrations are tabulated in the first column of table 5.4. At t=5 seconds, the 
HHC loops are engaged and the nominal T-matrix controller begins to reduce the 4/rev 
vibrations to a lower level. It takes approximately 2-3 seconds for the 4/rev vibrations 
to reach a new steady state condition where 67% of 4/rev in-plane vibrations have been 
reduced (Table 5.4). The large time constant of 2-3 seconds consists with the slow HHC 
loop dynamics stated in section 5.2.1. 

5.2.4 Transient vibration in maneuvering flight 

While the impact of HHC on handling-qualities is negligible, there are significant vibration 
responses to pilot inputs in both the baseline (HHC-off) case and the nominal (k= 1) 
case. Figures 5.15 and 5.16 show the large transient responses for a -50° roll maneuver 
(moderate) starting from t=12 seconds. Once the maneuver is completed, the vehicle 
reaches a new trim vibration level. Similar results can also be observed in figures 5.17 
and 5.18, which demonstrate the large transient responses for a 20° pitch maneuver starting 
from t=12 seconds. 

Using F Xac as an example, figure 5.19 shows the A F XiC response of both the baseline 
and nominal case for the same -50° roll maneuver. The symbol A denotes the steady state 
vibration of F XiC at t=12 seconds has been removed from the figure. With the baseline case, 
figure 5.19 shows that there is a maximum transient peak excitation of 150 lb above the 
steady state vibration level in the F x iC channel. Note that the f x 4C steady state vibration 
level for the baseline case is 151 lb (Table 5.4). Therefore, this maximum transient peak 
excitation is roughly the same as the baseline steady state vibration level. With the nominal 
T -matrix controller engaged, the maximum F x 4C vibration transient increases to 163 lb, 
which is 9% higher than the baseline case. In other words, with the nominal T-matrix 
controller engaged, the transient vibration response during maneuvering flight reaches 
similar levels to the trim condition with HHC-off. Nevertheless, the nominal T-matrix 
controller is able to reduce the transient load back to lower levels faster than baseline case 
after the 15-second point. 

The performance of the HHC system in suppressing the vibration response to pilot input 
is also reflected in the frequency-responses: F X4C /8i at , F XiS /Si a t, Fy AC /5i a t, etc. The RMS, 
determined from the integral under the frequency-response squared functions, is a useful 
measure of the vibration response to the broadband pilot inputs for different HHC system 
designs. The spectral integration to determine the RMS is conducted up to a frequency of 
3 rad/sec. The 3 rad/sec cut-off frequency corresponds to the roll command bandwidth, 
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and it is a good estimate of the maximum closed-loop piloting frequency. Finally, the RMS 
levels were normalized using the baseline vibration RMS for the roll maneuver to show the 
relative improvement (or degradation) in vibration suppression by the HHC system. 

Figure 5.20 shows the frequency response of f x 4C with respect to the lateral pilot input. 
Looking at the magnitude curve, the nominal case has a small magnification effect at higher 
frequency range (1-3 rad/sec) and a large reduction effect below the frequency of 0.9 
rad/sec. Both effects are consistent with the result shown in figure 5.19, where there is 
a small increase in the transient vibration excitation and a large reduction in steady state 
vibration level. Because the nominal T-matrix controller is capable of suppressing the 
F X/IC / dint vibration response more than it magnifies, there is 4.1% reduction in F X4C /5i at 
channel. The small reduction of 4.1% does not seem to reflect amount of vibration 
suppression shown in the figure. This is because the figure is on the logarithm scale, 
which biases toward the lower range. When including other seven channels ( Fx 4s /Si a t , 
F Y4C /Siat, Fy 4S /6i at , Fx 4C /5i on , Fx 4S /Si on , F Y4C /$i on , F Y4S /5 lon ), the average vibration in 
maneuvering flight for a nominal case is 3.2% above the baseline case (Table 5.5). This 
shows that the nominal T -matrix controller is ineffective for vibration suppression during 
maneuvers. 

5.2.5 Ideal integrator approximation 

Many previous studies (refs. 6, 18-20,22,23,61-64) represented the helicopter plant model 
in figure 5. 13 by a fixed T -matrix, which is a linear approximation of the vibration response 
to the HHC inputs at a steady-state condition. In other words, T-matrix corresponds to the 
linear state-space model at DC gain 1 to within the accuracy of the linear model extraction 
process. This method eliminates the need for a detailed model of the periodic helicopter 
dynamics. The nominal (k= 1) T-matrix controller (HHC Controller in fig. 5.13) is simply 
a k/s diagonal compensator multiplied by the fixed-gain regulator T^. The broken-loop 
response matrix (k/s) T t T will thus be a nearly diagonal matrix of k/s responses. 
This corresponds to single-input/single-output loop and without loop interactions (e.g., no 
response of the 3s loop to 3c transients). Assuming a nominal gain of k= 1, this ideal 
approximation gives loop crossover frequencies of u; c = 1 rad/sec, 90° phase margin, and 
infinite gain margin in every loop as illustrated in figure 5.21. 

Next, the helicopter vibration model is replaced with the LTI-HHC model. The actual 
broken-loop response for the 3/rev-cosine loop shown in figure 5.21 confirms that the k/s 
approximation is accurate for frequencies up to about the 1 rad/sec crossover frequency. 
There is a gain offset associated with the deviation between the steady response of the 
nonlinear simulation (T-matrix) and the steady- state response of the linearized model. For 
frequencies above 1 rad/sec, there is significant deviation from the 1/s ideal response, 
especially in phase, due to the dynamics of the 4/rev vibration response relative to the 
simple steady-state approximation (T-matrix). 


'DC Gain is the ratio of the output /input signal at the steady-state condition. 
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The Fx 4C vibration response to a unit pulse input is shown in figure 5.22 to be well 
damped. Increasing the HHC feedback gain (k= 2) raises the broken-loop crossover 
frequency and the closed-loop HHC disturbance rejection bandwidth (fig. 5.23). There 
is an associated reduction in the closed-loop transient settling time, as was also concluded 
by Shin et al. (ref. 5). But, there is also a magnification of peak disturbance at frequencies 
above crossover (fig. 5.23), which is consistent with classical control theory and which 
shows up in the time-domain as well (fig. 5.22). 

5.2.6 Optimized HHC controller 

Analyses with CONDUIT® show that an improvement in the suppression of vibration 
transients during the maneuvering flight can only be achieved by increasing the HHC 
crossover frequency to a value that is close to the 3 rad/sec piloted bandwidth. At 
this increased crossover frequency, the use of the T-matrix (which is a steady-state 
approximation) to simulate the helicopter vibration model is unacceptable for controller 
optimization and analysis, and must be replaced with the complete dynamic LTI-HHC 
model developed in chapter 4. Furthermore, the simple k/s HHC controller architecture 
must be augmented with the addition of a second order lead-lag compensator (fig. 5.24) in 
each loop to add robustness and achieve the needed stability margins. The HHC feedback 
controller now takes the form: 


Each HHC control loop contains five design parameters, and the same controller is used for 
the cosine and sine loops of a particular harmonic. Thus, for the three harmonics (6 loops), 
there are 15 HHC feedback parameters in total. 

CONDUIT® was used for HHC controller analyses and optimization. The key HHC 
design specifications included in the analysis were HHC loop stability margins and 
vibration suppression performance. Gain and phase stability margins were determined 
for each of the six broken HHC loops, and the vibration suppression performance are 
determined from the RMS value. The design metrics are list below (Appendix A): 

• Eigenvalue real part (EigLcGl) 

• Stability margins (StbMgGl) 

• Actuator RMS value (RisAcGl) 

CONDUIT® quickly minimized the sum of the normalized vibration RMS values for the 
four in-plane shears to both lateral and longitudinal input without sacrificing the required 
HHC loop stability margins. The optimum HHC feedback parameters are presented in 
table 5.6. The final evaluations of these HHC design specifications are shown in figure 5.25. 
Subfigures 5.25e, f, g, and h show the relative improvement/degradation in vibration 
suppression by the optimized HHC controller. The RMS> 1 represents the vibration level as 



(5.1) 
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increased with respect to the baseline (HHC-off) case, the RMS<1 indicates the vibration 
level as reduced with respect to the baseline case, and the RMS=1 represents the vibration 
level as the same as the baseline case. Except for the F Yac / Si xm channel, the vibration levels 
of other channels have been reduced. 

Following the previous example, figure 5.26 shows the frequency response of Fx 4C with 
respect to the lateral piloted input. The magnitude plot (top figure) shows that the optimized 
HHC controller has dramatically reduced the vibration response by 64% over broadband 
pilot lateral inputs. In terms of overall performance, the average vibration in maneuvering 
flight for the optimized HHC controller is 37% below the baseline case (Table 5.5). This 
is achieved by increasing the crossover frequencies to their maximum values (e.g., tu c = 
2.5 rad/sec in the 3/rev-cosine loop, Table 5.7) while still maintaining adequate stability 
margins (fig. 5.27). 

Similar conclusion can also be drawn from the time domain results. Figures 5.28-5.31 
are the time history of the vibration responses with the optimized HHC controller. The 
vibration responses of the nominal (fc=l) and baseline cases (HHC-off) are also presented 
in the figures. Fooking at F Xac / di nt in figure 5.28, the vehicle starts from a steady state 
condition, and maintains at a steady 4/rev vibration level. At t=5 seconds, the HHC loops 
are engaged and the optimized HHC controller begins to reduce the 4/rev vibrations to a 
lower level. Although the optimized HHC controller has reached the same new steady- 
state condition as the nominal T-matrix controller, the optimized HHC controller has a 
much lower raise time which is directly related to the higher crossover frequency. The 
peak vibration in Fx 4C /Slat channel shown in figure 5.32 is now 73 lb, or 51% below the 
baseline result, which again tracks the frequency-domain results of table 5.5 closely. One 
can clearly see that the optimized controller has achieved performance superior to that of 
the baseline (HHC-off) and nominal T-matrix controller cases. 
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Table 5.1. HHC broken-loop stability margins; nominal T-matrix controller. 


Broken-Loop 

Channel 

LU C 

(rad/sec) 

Gain Margin 
(dB) 

Phase Margin 
(deg) 

3/rev COS 

0.92 

16.0 

74.3 

3/rev SIN 

1.01 

14.9 

75.4 

4/rev COS 

1.03 

15.9 

75.1 

4/rev SIN 

0.99 

16.2 

75.9 

5/rev COS 

1.00 

15.8 

76.0 

5/rev SIN 

0.93 

16.5 

76.7 


Table 5.2. Second order actuator model parameters. 


Nature Frequency, u, (rad/sec) 

30.0 

Damping Ratio, ( 

0.8 

Rate Saturation Limit (in/sec) 

600.0 

Upper Position Limit (in) 

60.0 

Lower Position Limit (in) 

-60.0 


Table 5.3. Flight control system parameters. 


K u 

0.000 

K u 

0.000 

K v 

0.000 

K p 

-0.570 

K q 

-0.766 

K r 

-0.693 

k 4> 

-2.480 

K e 

-3.416 

K 4> 

-2.565 

KI,p 

1.084 

KIg 

0.000 

KI W 

0.000 
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Table 5.4. Effect of fixed T-matrix on steady state vibration level. 



Baseline 

(HHC-off) 

Nominal (fc=l) 
T -matrix Controller 

Percent 

Changed 

Fx ac (lb) 

151.6 

51.4 

-66.1% 

F XiS (lb) 

87.8 

21.7 

-75.3% 

Fy 4C (lb) 

73.5 

-3.4 

-95.4% 

Fy as (lb) 

-61.3 

-42.6 

-30.5% 


Average -66.8% 


Table 5.5. Vibration RMS with respect to piloted roll and pitch inputs. 


Nominal T -matrix Optimized HHC 
(k= 1) (Lead-Lag) 


Roll Maneuvering Llight 


Fx 4C 

-4.1% 

-63.8% 

F X 4S 

1.4% 

-21.3% 

Fy 4C 

20.1% 

-15.4% 

f y 4S 

13.6% 

-67.9% 

Pitch Maneuvering Plight 

F XiC 

-1.8% 

-46.4% 

Fx 4S 

-59.9% 

-51.1% 

f y 4C 

43.5% 

9.7% 

F y 4S 

13.0% 

-40.8% 

Average 

3.2% 

-37.1% 


Normalized relative to the baseline RMS 
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Table 5.6. HHC controller parameters. 



Nominal T -matrix 

Optimized Lead-Lag 


Controller 

Controller 

K 3p 

1.000 

1.153 

K± p 

1.000 

1.663 

K-bp 

1.000 

1.236 

^ 77.3 


1.463 

UJnA 


5.253 



2.848 

Cn3 


2.539 

Cn4 


0.583 

Cn5 


1.246 

Ud3 


6.900 

^dA 


5.627 

Ud5 


6.787 

Cd3 


3.494 

CdA 


1.207 

Cd5 


1.179 


Table 5.7. HHC broken-loop stability margins; optimized HHC controller. 


Broken-Loop 

Channel 

LU C 

(rad/sec) 

Gain Margin 
(dB) 

Phase Margin 
(deg) 

3/rev COS 

2.46 

8.2 

74.6 

3/rev SIN 

2.71 

6.0 

78.6 

4/rev COS 

1.52 

17.6 

52.8 

4/rev SIN 

1.45 

18.0 

55.5 

5/rev COS 

1.80 

6.0 

99.6 

5/rev SIN 

1.33 

6.6 

102.3 
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Phase (deg) Magnitude (dB) 




Figure 5.1. Analytic model validation, baseline (HHC-off) case. 
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Phase (deg) Magnitude (dB) 




Figure 5.2. Effect of HHC input on rigid body dynamics, nonlinear model. 
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Phase (deg) Magnitude (dB) 




Figure 5.3. Effect of HHC input on rigid body dynamics, LTI-HHC model. 
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Figure 5.4. General closed-loop HHC vibration reduction scheme. 
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(b) LTI-HHC model (Continuous time-domain system) 


Figure 5.5. Closed-loop HHC systems comparison. 




(a) 0 3C broken-loop responses comparison. 




(b) 0 3S broken-loop responses comparison. 


Figure 5.6. HHC 3P broken-loop frequency responses comparison. 
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Phase (deg) Mag (dB) Phase (deg) Mag (dB) 




(a) 9 ac Broken-loop responses comparison. 




(b) 9 as Broken-loop responses comparison. 


Figure 5.7. HHC 4P broken-loop frequency responses comparison. 
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(a) 0 5C Broken-loop responses comparison. 




(b) 0 5S Broken-loop responses comparison. 


Figure 5.8. HHC 5P broken-loop frequency responses comparison. 
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Figure 5.9. HHC vibration reduction system, HHC-loops disengaged. 
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#3 Actuator 
command 

-► 

#4 Actuator 
command 


Figure 5.10. Automatic flight control system schematics. 
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Actuator 

Response 


Figure 5.11. Second order actuator schematics. 
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(b) EigLcGI: 
Eigenvalues (All) 
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Figure 5.12. CONDUIT® handling-quality design specifications; HHC-loops disengaged. 
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Figure 5.13. HHC vibration reduction system, HHC-loops engaged (k= 1). 
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Figure 5.14. CONDUIT® handling-quality design specifications; HHC-loops engaged 
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Figure 5.15. F\ vibration response in roll maneuvering flight, T-matrix controller, nominal 
case (k= 1). 
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Figure 5.16. Fy vibration response in roll maneuvering flight, T -matrix controller, nominal 
case (k= 1). 
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Figure 5.17. Fx vibration response in pitch maneuvering flight, T-matrix controller, 
nominal case (k= 1). 
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Figure 5.18. Fy vibration response in pitch maneuvering flight, T-matrix controller, 
nominal case (k= 1). 
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Figure 5.19. f x 4C vibration response in roll maneuvering flight, T-matrix controller, 
nominal case (k= 1). 


146 


Phase (deg) Magnitude (dB) 




Figure 5.20. Fx AC vibration response in roll maneuvering flight, T-matrix controller, 
nominal case (k= 1). 
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Figure 5.21. Ideal and actual broken-loop responses comparison, 3/rev-cosine loop, 
nominal T-matrix controller case. 
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Figure 5.22. F Xac unit pulse response; T-matrix controller, k= 1 and 2. 
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Figure 5.23. Broken-loop response of 3/rev input; T-matrix controller, /;:= I and 2. 
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Figure 5.24. HHC vibration reduction system, HHC-loops engaged, optimized lead-lag 
compensator. 
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Figure 5.25. CONDUIT® HHC design specifications; HHC-loops engaged, optimized 
lead-lag compensator. 
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Figure 5.26. f x ac vibration response in roll maneuvering flight, optimized lead-lag 
compensator. 
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Figure 5.27. Broken-loop response of 3/rev input; optimized lead-lag compensator. 
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Figure 5.28. Fx vibration response in roll maneuvering flight; optimized lead-lag 
compensator. 


155 


otA d (qi) st,A d 




Figure 5.29. Fy vibration response in roll maneuvering flight; optimized lead-lag 
compensator. 
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Figure 5.30. Fx vibration response in pitch maneuvering flight; optimized lead-lag 
compensator. 
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Figure 5.31. Fy vibration response in pitch maneuvering flight; optimized lead-lag 
compensator. 
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Figure 5.32. F XiC vibration response in roll maneuvering flight; optimized lead-lag 
compensator. 
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6 Summary and Conclusions 

The increasing opportunities provided by novel sensing and actuation technologies, and 
the advancements in the theory and practice of flight and rotor control systems open 
unprecedented possibilities in the constant search for low vibration levels and favorable 
handling qualities in modem helicopters. At the same time, however, greater care than ever 
must be taken to ensure that these advanced controls cooperate harmoniously and prevent 
adverse dynamic interactions. 

The present work makes a contribution toward this goal by developing new mathematical 
tools for the analysis and design of active rotor control systems, more specifically, Higher 
Harmonic Control (HHC) systems, and by using these tools to carry out the first systematic 
study of the interaction of HHC and Automatic Flight Control Systems (AFCS) available 
in the literature. 

This chapter provides a summary of the work presented in the research, details 
conclusions drawn from its results, and outlines some recommendations for future work. 
Chapter 2 describes the key features of the formulation and solution techniques for 
the baseline helicopter simulation model used in this study. Chapter 3 provides basic 
information on the HHC algorithm. The extraction of a linearized, time-invariant dynamic 
model of the helicopter that includes higher harmonic content is a key contribution of this 
study, and is described in detail in chapter 4. Other important contribution, namely, the 
AFCS -HHC interaction study, is presented in chapter 5. AFCS design procedures, and 
basic concepts of Fourier analysis and treatment of rotor degrees of freedom, are briefly 
reviewed in the Appendices. 

6.1 Summary 

A realistic analysis of the interaction between AFCS and HHC requires a mathematical 
model of a helicopter of adequate sophistication. This model must be able to provide 
sufficiently accurate predictions of vibratory loads in both trimmed and maneuvering 
flight. This model was described in chapter 2. An existing, state-of-the-art flight dynamic 
simulation model was improved to allow the calculation of vibration levels both at the 
center of mass of the helicopter and at specific locations such as pilot and copilot seats. 
The results obtained with this model were successfully validated through comparison with 
other simulation models and with flight test data. 

An HHC system is composed of several elements which must all be modeled in a 
rigorous mathematical way. This is the main topic of chapter 3. The harmonic analyzer, 
which extracts the desired frequency components of the rotor vibrations, was studied first. 
A Fourier analysis method was described, and the effects of windows were discussed. 
Then, the HHC control algorithm was presented, in the traditional T-matrix form, and 
validated through simulation. Finally, issues associated with the discrete, rather than 
continuous, implementation of HHC were discussed. 
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The methodology for the extraction of a high-order, time-invariant linearized model of 
the coupled rotor-fuselage system was systematically described in chapter 4. This model 
included both the pilot and the HHC inputs, and both the averaged and the high frequency 
dynamics of the rotor states. The resulting model contains, as a subset, the more traditional 
linear time-invariant representation without high frequency rotor dynamics and higher 
harmonic controls. Therefore, the description of the methodology started with this well- 
known subset. The methodology to extract the remaining partitions of the state, control, 
and output matrices was presented next, with partitions chosen in an order to allow the 
progressive introduction of the new key concepts. Chapter 4 continued with the application 
of a technique to a simplified rotor model, entirely formulated analytically. This model 
was not sufficiently sophisticated to be used in the remainder of the research. However, it 
was very useful to illustrate and validate the methodology. In fact, the higher harmonics 
of the rotor motion and of the control inputs were explicitly accessible in the equations in 
analytic form. A more complete validation, performed by comparing hub loads and rotor 
states predicted by the linearized model and by the full nonlinear simulation, concluded the 
chapter. 

The newly developed linearized model was then used to carry out a study of the 
interaction between HHC and AFCS, described in chapter 5. First, the effect of open- 
loop HHC on rigid body dynamics was examined in detail, by observing the changes in the 
frequency responses of helicopter to pilot inputs when the HHC controller was turned on. 
Then, a full closed-loop interaction study was performed. The study included a validation 
through simulation of the response of the helicopter with all control loops closed, an 
analysis of the vibratory loads with and without HHC in both trimmed and maneuvering 
flight, and a discussion of the tailoring of the HHC controller to improve its performance 
in transient maneuvers. 


6.2 Conclusions 

This section presents the main observations originating from this research, and the key 
conclusions of the study. The conclusions related to the new linearization procedure are 
presented first, followed by those concerning the AFCS/HHC interaction study. 


6.2.1 Extraction of linearized, time-invariant models 

1. The traditional constant coefficient linearized models of coupled rotor-fuselage 
dynamics, obtained through multiblade coordinate transformations followed by 
averaging over one rotor revolution, are not suitable for studies involving rotor 
vibrations, even if the control vector includes the higher harmonics typical of HHC. 
In fact, the averaging removes all higher harmonics of the rotor response. Such a 
model will capture the effects of HHC on the low frequency rigid body motion of the 
helicopter and of the tip path plane, but not on the N/rev vibrations. 
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2. The constant coefficient linearized model developed in this research, which explicitly 
includes states describing the high frequency rotor dynamics, does capture the 
vibratory loads, and the effects that HHC can have on them. The price for modeling 
vibrations with a linear time-invariant system, compared with a linear system with 
periodic coefficients, is an increase in the size of the system. On the other hand, the 
entire arsenal of tools of linear time-invariant system theory can now be used. 

3. The validation with the full nonlinear simulation model shows that there is very good 
agreement between the hub loads predicted by the new LTI-HHC model and the hub 
loads of the nonlinear model, both for HHC and pilot inputs. This suggests that 
a linearized model that intrinsically includes higher rotor harmonics is sufficiently 
accurate for full load predictions, at least for the aircraft configuration and flight 
conditions considered in this study. In other words, periodicity plays a far more 
important role than nonlinearity. 

4. One limitation of the current LTI-HHC model is that it can model only the 4/rev 
components of the system and not the higher frequency components that enter the 
fuselage, i.e., 8/, 12/rev, etc., for the 4-bladed rotor of this study. However, this 
limitation can be easily overcome, by including additional harmonics in the LTI- 
HHC model using the same methodology as for the 4/rev states. 

5. Possibly for historical reasons, the starting point for the vast majority of HHC 
modeling research and applications has been an update equation that links the 
vibration harmonics to the HHC harmonics through the T-matrix. Using instead 
an ( A , B, C, D ) state-space representation as a starting point leads to a much richer 
and informative picture. In fact, the traditional update equation is included as a 
subset (through a partition of the control matrix B ), and the additional effects on 
vibrations of pilot inputs and all the states, including aircraft rigid body, rotor, and 
inflow states, are now modeled explicitly. These additional effects are not included in 
the traditional update equations, and are usually taken into account indirectly through 
on-line identification and adaptation schemes. 

6.2.2 HHC/AFCS interaction study 

1 . In general, the closed-loop HHC system has little influence on the handling qualities 
characteristics of the helicopter, and on the behavior of the flight control system, at 
least for the articulated rotor configuration used in this study. This conclusion is 
drawn on the basis of the analysis of the effects of HHC on the frequency response 
to pilot inputs. The effects of HHC on trim were not addressed explicitly, but 
the simulated free flight responses with HHC suggest that these effects are not 
significant. 

2. Although the typical 3/, 4/, and 5/rev HHC inputs for a 4-bladed rotor are at high 
frequency (81, 108, and 135 rad/sec, respectively, for the helicopter used in this 
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study), the crossover frequency of each HHC loop is only about 1 rad/sec. Because 
of the high HHC input frequency, one might expect a large frequency separation 
between the flight control and the HHC inputs, and assume that these two systems 
would not interfere with each other. Instead, the results clearly show that this is not 
the case, and that the potential for AFCS/HHC interaction does exist. 

3. The vibration response to maneuver inputs, and not just to steady state inputs, must 
be considered as part of the HHC system design process. If the HHC algorithm is not 
properly designed, the transient vibrations in the early phases of a maneuver could 
be higher than if no HHC system was present. 

4. An HHC controller that improves the suppression of vibration transients has the 
higher loop crossover frequencies. For the cases studied, these frequencies are of 
the order of 3 rad/sec. At these frequencies, the use of the T-matrix approach to 
simulate the helicopter vibration model is unacceptable for controller analysis and 
optimization. This is because the T-matrix is simply a k/s diagonal compensator 
multiplied by a fixed-gain regulator, and a comparison with the more sophisticated 
LTI-HHC model developed in this study shows that it is inaccurate for crossover 
frequencies greater than about 1 rad/sec. Increasing the T-matrix controller 
feedback gain (k= 2) reduces the closed-loop transient settling time and increases 
the magnitude of the peak disturbance at frequencies above crossover frequency. 

5. For the maneuvering flight conditions considered in this study, the optimized HHC 
system designed using the new linearized model reduces vibratory hub shears by 37% 
compared to the baseline case, and 39% compared to nominal T-matrix controller 
case. Therefore, the need for on-line identification and adaptation of the T-matrix is 
greatly reduced if not completely eliminated. This is important from a practical point 
of view, because of the danger that an adaptive system on board a helicopter might 
react in unpredictable and unwanted ways, which can clearly create safety-of-flight 
issues. 

6.3 Future work 

The research presented in this study has shown the importance of the HHC/AFCS 
interaction on the transient vibration suppression. However, there are areas in which the 
present analysis was limited. This section suggests some areas for improvement. 

1. Improve the flexible blade model, for example by adding additional blade modes 
and increasing the number of blade finite elements. While the vibration results can 
be considered qualitatively representative, a more sophisticated model is probably 
needed for quantitative evaluations (e.g., for a precise quantification of the benefits 
of HHC). 
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2. For the same reason just mentioned, improve airload calculations, especially by 
adding non-uniform inflow and unsteady aerodynamics modeling. For the design 
of the HHC system, the improved models must obviously be in state-space form (not 
necessarily linear). This characteristic is not required for validation purposes. 

3. Further validate the vibratory hub load level predicted by the mathematical model 
with wind tunnel data or flight test data. Because of the large scatter, the flight test 
data used in this study were only adequate for a qualitative validation. Unfortunately, 
no other flight test data was publicly available for a helicopter without some or all of 
the normal vibration suppression devices. 

4. Repeat the study with a helicopter configuration with lowly damped coupled 
rotor/body modes, such as hingeless or bearingless rotor helicopters. The articulated 
rotor configuration used in this study had hydraulic lag dampers, and aeromechanic 
stability was never an issue. 

5. Apply advanced control design theories such as H 2 , Hoc control design methods in 
an attempt to achieve further improvements in vibration reduction that may remove 
the need for adaptive T-matrices. 
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APPENDIX A 


CONDUIT HQ-WINDOW SPECIFICATIONS 


Following are the handling-qualities specifications and control system metrics used in the 
AFCS optimization procedure: 

(a) Eigenvalues (All) This criterion is used to ensure that all the real parts of the 
eigenvalues of the system are zero or negative, ensuring that all the dynamics are stable 
or neutrally stable. At any given iteration, the sum of unstable eigenvalues real parts or the 
largest stable eigenvalue is returned as the spec metric. 

(b) Minimum Crossover Frequency The crossover frequency is defined as the 
frequency where the magnitude curve crosses 0 dB. For multiple crossings, the highest 
crossover frequency is returned. This specification is intended as a hard constraint to a 
greater than zero value of crossover frequency. 

(c) Gain/Phase Margins The spec has very sophisticated logic for treating stable, 
conditionally stable, and unstable systems. It also has logic for correctly accounting for 
right-half plane poles and zeros. A table of margins is built for all crossings of the 0 dB and 
-180 deg lines and displayed in the supporting plot. The spec returns the minimum gain and 
phase margin values from the table. The level 1 boundaries are taken from MIL-F-9490D. 

(d) Bandwidth Specification The vehicle response to cockpit control force or position 
inputs shall meet the limits specified. It is desirable to meet this criterion for both controller 
force and position inputs. If the bandwidth for force inputs falls outside the specified limits, 
flight testing should be conducted to determine that the force feel system is not excessively 
sluggish. 

(e) Attitude Response Damping Ratio (from peak overshoot) The calculation of the 
damping ratio (zeta) is from peak overshoot of the time response to a step input. ADS-33D 
required a minimum damping ratio of 0.35. Systems whose eigenvalues all have damping 
ratios of greater than 0.35 could still have excessive overshoot due to the presence of zeros 
in the response. This spec ensures that the end-to-end attitude response has an effective 
damping ratio greater than 0.35 base on the time response. An appropriate input should be 
used to result in a step response. 

(f) Crossover Frequency The crossover frequency is defined as the frequency where the 
magnitude curve crosses 0 dB. For multiple crossings, the highest crossover frequency is 


167 


returned. This specification is intended as an objective to minimize crossover frequency in 
CONDUIT® phase 3 optimization. 

(g) Damping ratio This specification is used to ensure damping is above the minimum 
value specified. This is achieved by checking the damping ratios of the eigenvalues within 
the range of natural frequencies specified. 

(h) Rise Time (Calculated from 10% to 90% of peak response) This spec estimates 
rise time for first-order SISO systems by finding the peak of the time response, and 
calculating the time between 10% and 90% of the peak magnitude. 
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APPENDIX B 


FOURIER TRANSFORMS 
B.l Fourier Transform (FT) 

Fourier transform can be viewed as a generalization of the Fourier series representation of a 
periodic function. Unlike the Fourier series which is an approximation of the source signal, 
Fourier transform is a direct mapping between time-domain and the frequency-domain, and 
it is fully reversible. 

Let /(f) be a continuous-time signal, its continuous Fourier transform F(u) is defined 
by 

/ OO 

f(t) e~^ ut dt, — oo < to < oo (B.l) 

- OO 

where u is the frequency variable in rad/sec. In many applications, the source signal 
/(f) cannot be given in common function 1 ; therefore, Fourier transform is often computed 
numerically. This numerical computation can be performed in either the continuous-time 
domain (continuous Fourier Transform) or the discrete-time domain (discrete-time Fourier 
transform). 

Because a digital computer works only with discrete data, numerical computation of the 
Fourier transform of /(f) requires discrete sample values of /(f). In addition, a digital 
computer can compute the transform F(u) only at discrete values of cu; therefore, discrete- 
time Fourier transform is often used in many applications. 

B.2 Discrete-Time Fourier Transform (DTFT) 

Let / (k) be a sampled version of a continuous-time signal / (f) with f evaluated at sample 
time f = fcT, where T is the sample interval. 

f(k) = f(t) \ t = kT = /(AT), k = 0, ±1, ±2, ... (B.2) 

The Fourier transform of / (k) is defined by 

OO 

= £ f(k)e — oo < < oo (B.3) 

k=— oo 

Note that DTFT is directly analogous to the FT, and it is not an approximation to the FT. 

The DTFT requires the calculation of the sums of equation B.3 for all frequencies range. 
In practice, F (Q.) is usually computed only for a discrete set of frequency variable O, and 


'it is the generalized transform typically shown in Fourier transform table 
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this is accomplished by using the N-point discrete Fourier transform (N-point DFT). 

N - 1 

F n = f (k) e- j2 " fcn/N , n = 0, 1,..., N - 1 (B.4) 

k = 0 

where N is a positive integer. 

B.3 Fast Fourier Transform (FFT) 

The computation of equation B.4 can be carried out using a fast algorithm called the 
Fast Fourier Transforms. It is a new N-point DFT algorithm developed by Tukey and 
Cooley (ref. 65) in 1965 which reduces the number of computations from something on the 
order of N 2 to N logiV. Because many special computers or add-on cards are available to 
perform the FFT algorithm at ultra-high speed, FFT opens the possibility of a wider use 
of the FT in many other areas such as the computational physics and many engineering 
applications. Additional information regarding to FT, DTFT, DFT, and FFT can be found 
in references 66 and 67. 
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APPENDIX C 


MULTI-BLADE COORDINATE TRANSFORMATION 


C.l Converting rotor equations of motion 

Let the blade flapping equations of motion for a four-bladed rotor in the rotating system be 

xr + C R ± R + K R ~x. R = i R (C.l) 


where 


Xfi 


1 T 


Pi i 02 1 03 i Pa 


(C.2) 


The matrix Tp is the multi -blade coordinate transformation which converts x from the 
fixed to the rotating system as follows: 


rji R 

Xfl = T F Xp 


(C.3) 


The first and the second time derivative of equation C.3 are: 


x B = TpX F + Tp± F (C.4) 

xr = Tp xp + 2 Tp itp + Tp xp (C.5) 


Substitute equations C.4 and C.5 into equation C.l yields 

Tp Xr + 2Tp Xr + T# Xr + Cr (T# Xr + T# Xr) + Kr T# Xr = f* (C.6) 
Multiply Tp through equation C.6 and re-arrange the equation, equation C.6 becomes 


where 


Xp + CrXr + ArXr = fR (C.7) 


1 T 


Xr 

— Po, Pic , /?ls, P 2 

(C.8) 

Cr 

= Tp (C R T* + 2 T#) 

(C.9) 

A> 

= Tp (TP + C R fP + Ar Tj?) 

(C.10) 

fR 

= U^R 

(C.ll) 
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(C.14) 
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(C.16) 


(C.17) 


C.2 Converting state-space representation 

Let equation C.18 be the state-space representation of rotor equations of motion in rotating 
system. 

± R = A r x r + B r u (C.18) 
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where A R and B R is the state matrix and the control matrix in the rotating system, 
respectively. Substituting equations C.3 and C.4 into equation C.18 yields: 


Tp Xp + Tp yip — A r Tp -x.p + B r u 

(C.19) 

Multiply Tp through equation C.19 and re-arrange the equation as: 


ip = Ap ^.p + Bp u 

(C.20) 

where 


= t r (a r t# - tp) 

(C.21) 

Bp = Tp B r 

(C.22) 
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